Calculando el Spread y el Bounty descritos por Brynjolfsson & McAfee mediante PCA y Factor Analysis
INTRODUCTION
El objetivo de este trabajo es realizar un análisis multivariante que nos permita indagar en la naturaleza de dos variables descritas por Brynjolfsson & McAfee (2014) en The Second Machine Age, como son el spread y el bounty. En concreto el spread creemos que será más sencillo de encontrar mediante un anáisis factorial que el bounty, ya que este último tiene una descripción algo más confusa. Podemos utilizar la descripción de esta sinopsis para definir ambas, pues resume muy bien ambos conceptos:
‘What the authors call “bounty” is their attempt to measure the benefits of new technology in ways reaching beyond such measures as GDP, which they say is inadequate. They use “spread” as a shorthand way to describe the increasing inequality that is also resulting from widespread new technology.’
Además, una vez hayamos conseguido definir bien la variable latente spread, podremos ver cómo están posicionados los países del mundo respecto a la misma, y cuáles están desmarcándose del resto positiva o negativamente debido al progreso tecnológico y a una creciente desigualdad en la repartición de los beneficios del mismo.
Para ello utilizaremos datos de distintos datasets, entre los cuales están UN, WHO, WB. En todos estos datasets existe una gran ausencia de información, una infinidad de valores ausentes y apenas son capaces de ponerse de acuerdo para utilizar los mismos nombres para los países. Todo esto dificultará la obtención del factor óptimo que nos de el spread (o la posición dentro de esa escala de spread de cada uno de los países).
Para tratar de hallar este spread, se utilizarán variables económicas (las que se pueda, ya que muchas que además son interesantes como el GINI ratio tienen muchos valores ausentes – una pena porque esta variable apunta en la misma dirección de los autores, hacia la desigualdad, nos da una idea de cómo es la desigualdad dentro de un mismo pais, lo cual nos podría ayudar para sacar tanto el spread como el bounty), de salud (acceso a sanidad y al agua, mortalidad infantil…), de igualdad (porcentaje de mujeres en el labour force…), etc. Buscamos tener variables que no midan sólo el potencial económico de un país, ya que como bien dicen los autores este potencial económico medido en términos de GDP es insuficiente a la hora de explicar el progresivo spread que se está generando entre unos países y otros en esta Second Machine Age. Habría sido muy positivo incluir variables que reflejen progreso tecnológico, venta de patentes etc; de hecho estas variables fueron descargadas del World Bank, pero tienen una infinidad de valores ausentes por lo que serán eliminadas al menos en su mayoría como puede verse en el código más adelante.
Otras medidas generales de progreso que investigaremos serán lo ecológico que es cada país en el uso de las energías, ya que como bien explican los autores, este “spread” debido al progreso tecnológico tendrá que ser sostenible para poder considerarse como tal, pues en el largo plazo no se mantendría de otra forma. El uso de internet y de celulares, la proporción de la economía dedicada a la agricultura, etc serán otras variables que se explorarán. Primero nos descargaremos todas las variables que nos encajen dentro de esta definición que dan Brynjolfsson & McAfee, bien sea para medir el spread o el bounty (aunque este ya decimos que es más difícil de medir, más aun sin poder contar con variables como el GINI ratio, que nos da la desigualdad de un país, que podría entenderse como la inversa de los beneficios del desarrollo tecnológico, tal y como lo describen ellos, ya que el bounty se encarga de explicar los beneficios para la población en general. De todas formas, creemos de partida que esta variable sólo podría entenderse como tal si pudiéramos hacer un análisis a lo largo del tiempo del factor que identifiquemos con la misma, ya que de otra forma no existe manera de constatar que efectivamente significa mejora de la calidad de vida debido al progreso tecnológico).
Después de descargar las variables y hacer un tratamiento exhaustivo de las mismas (por ejemplo armonizando los nombres y los códigos de los países en las diferentes bases de datos que nos bajemos), eliminaremos inteligentemente algunas observaciones; primero eliminamos las columnas que contienen demasiados valores ausentes, y después, sobre el data frame resultante, limpiamos los países que tienen información de pocas variables. De esta forma nos quedaremos con una matriz más completa, con un número inferior de valores ausentes, por lo que al utilizar la librería mice para la imputación de los mismos el random forest que utilizaremos para imputar variables tendrá suficientes observaciones cubiertas con respecto a las que tiene que rellenar. Al hacerlo de esta manera y no imputar directamente todo el dataset quedándonos con todas las variables que nos apetezca por muchos valores que tengan, reducimos la posibilidad de que se introduzca un sesgo en la muestra por haber utilizado pocas observaciones completas para rellenar muchos valores ausentes, en cuyo caso independientemente del modelo que uses te estás inventando los datos, con lo que el análisis completo quedaría invalidado.
Una vez tengamos el dataset completo con las variables y observaciones que queremos, procederemos a hacer un pequeño análisis descriptivo de las mismas; dibujaremos algunos gráficos, principalmente, para visualizar las distribuciones de variables, así como posibles outliers. También analizaremos la muestra de datos completa en su conjunto, viendo cómo se relacionan todas las variables con las demás. Por último, se realizará primero un PCA para ver cuánto porcentaje de la varianza podría coger con los primeros componentes, lo cual nos dará una idea también de qué tal seleccionadas están las variables. Posteriormente, trataremos de encontrar las variables spread y bounty realizando un Factor Analysis (toda la explicación respecto al número de factores que escoger etc se explicará más adelante). El modelo tiene la forma siguiente:
\(Variable_i = \mu_{i} + L_{i, j} \cdot Factor_j + \epsilon_{i}\)
\(i = 1, 2, ..., n\)
\(j = 1, 2, ..., m\)
n = número de observaciones
m = número de factor.
Finalmente se dibujarán algunos gráficos del mundo con los factores extraídos del análisis, y veremos si cuadra intuitivamente con la idea que tenemos de spread y bounty a partir de lo leído en The Second Machine Age
El trabajo está estructurado de la siguiente forma:
Recogida, tratamiento y selección de variables para el análisis.
Análisis descriptivo.
2.1. Análisis descriptivo individual.
2.2. Análisis descriptivo conjunto.
2.3. Inferencia.
- PCA & Factor Analysis: En la búsqueda del spread y el bounty.
3.1. PCA
3.2. FA
- Visualización de resultados,
1. RECOGIDA, TRATAMIENTO Y SELECCIÓN DE VARIABLES PARA EL ANÁLISIS.
En esta parte nos dedicaremos a leer y limpiar las variables. La mayoría de esta parte es código, ya que no hay mucho que comentar, pues no se centra en el análisis. En algunos casos cogeremos las variables para el último año, mientras que en otros casos tendremos que usar la media de los últimos años, o datos de 2016. Todo esto se verá explícitamente en el código. El procedimiento general es el siguiente: leemos la variable, utilizamos dplyr para dejarla en el mismo formato que a las demás, cambiamos nombres de columnas, seleccionamos solo las columnas que nos queremos quedar, identificamos de qué forma vienen los NAs para codificar valores como “…” como NAs, y hacemos merge para juntar cada variable a todas las demás.
Como hemos seleccionado sólo unas pocas de las variables totales consideradas, sólo las describiremos e indicaremos de dónde nos las bajamos de estas variables, para no alargar excesivamente sin sentido el trabajo.
Este dataset de aquí abajo está bajado de aquí, y contiene datos recopilados por las naciones unidas de todos los países del mundo. Nos quedaremos de ese dataset sólo con las variables que nos resulten de especial interés. Sin embargo, debemos decir que es un dataset mucho más completo que los encontrados en World Health Organization y en el World Bank; para las mismas variables tienen en la mayoría de los casos muchas más observaciones, no contienen tantos datos ausentes y es además un mapa más completo del mundo, al contener muchos más países.
Como los nombres de países no son homogéneos, trataremos de armonizar los nombres de estos en los distintos datasets de la siguiente forma:
Country_Name_correct<- c("Antigua and Barbuda"="Antigua",
"Bahamas, The"="Bahamas",
"Brunei Darussalam"="Brunei",
"Cabo Verde"="Cape Verde",
"Congo, Dem. Rep."="Democratic Republic of the Congo",
"Congo, Rep."="Republic of Congo",
"Cote d'Ivoire"="Ivory Coast",
"Egypt, Arab Rep."="Egypt",
"Faeroe Islands"="Faroe Islands",
"Gambia, The"="Gambia",
"Iran, Islamic Rep."="Iran",
"Korea, Dem. Rep."="North Korea",
"Korea, Rep."="South Korea",
"Korea, Republic"="South Korea",
"Republic of Korea"="South Korea",
"Kyrgyz Republic"="Kyrgyzstan",
"Lao PDR"="Laos",
"Macedonia, FYR"="Macedonia, Federal States of",
"Micronesia, Fed. Sts."="Micronesia, Federal States of",
"Russian Federation"="Russia",
"Slovak Republic"="Slovakia",
"St. Lucia"="Saint Lucia",
"St. Martin (French part)"="Saint Martin",
"St. Vincent and the Grenadines"="Saint Vincent",
"Saint Vincent and the Grenadines"="Saint Vincent",
"Syrian Arab Republic"="Syria",
"Trinidad and Tobago"="Trinidad",
"United Kingdom"="UK",
"United States"="USA",
"United States of America" = "USA",
"Venezuela, RB"="Venezuela",
"Virgin Islands (U.S.)"="Virgin Islands",
"United States Virgin Islands"="US Virgin Islands",
"Yemen, Rep."="Yemen",
"Sint Maarten (Dutch part)"="Sint Maarten",
"West Bank and Gaza"="Palestine",
"The former Yugoslav Republic of Macedonia" = "Macedonia",
"Democratic People's Republic of Korea"="North Korea",
"Bolivia (Plurinational State of)"="Bolivia",
"Lao People's Democratic Republic"="Laos",
"Micronesia (Federated States of)"="Micronesia, Federal States of",
"Venezuela (Bolivarian Republic of)"="Venezuela",
"State of Palestine" = "Palestine",
"China, Hong Kong SAR"="Hong Kong",
"China, Macao SAR" = "Macao SAR, China",
"Iran (Islamic Republic of)"="Iran",
"Korea, Dem. People's Rep."="North Korea",
"United Republic of Tanzania"="Tanzania",
"Republic of Moldova"="Moldova",
"Viet Nam" = "Vietnam",
"Virgin Islands"="US Virgin Islands",
"Hong Kong SAR"="Hong Kong",
"Micronesia" = "Micronesia, Federal States of",
"Czech Republic" = "Czechia",
"Macao" = "Macao SAR, China",
"Congo" = "Republic of Congo",
"Saint Kitts and Nevis"="St. Kitts and Nevis",
"Swaziland"="Eswatini")
correct_names <- function(df, correct_list = Country_Name_correct) {
df2 <- df
for (c in names(Country_Name_correct) ) {
df2[df2$country==c,"country"] <- Country_Name_correct[c]
}
return(df2)
}
add_code_absent <- function(df, countries = c("Eswatini", "Micronesia"), codes = c("SWZ", "FSM")) { #estos dos países nos dan fallos al meterles el código, así que usaremos la información de: https://www.iso.org/obp/ui/#iso:code:3166:SZ
for(i in 1:length(countries)) {
df[df$country == countries[i], "country_code"] <- codes[i]
}
return(df)
}
if(!require(countrycode)) {
install.packages("countrycode")
}## Loading required package: countrycodelibrary(countrycode)country_profile_vars <- read.csv("country_profile_variables.csv", sep = ",", fileEncoding = "utf-8", stringsAsFactors = F)
country_profile_vars <- correct_names(country_profile_vars)
country_profile_vars$country_code <- countrycode(country_profile_vars$country, "country.name", "iso3c")## Warning in countrycode(country_profile_vars$country, "country.name", "iso3c"): Some values were not matched unambiguously: Channel Islands, Eswatinicountry_profile_vars <- add_code_absent(country_profile_vars)
country_profile_vars <- country_profile_vars[!is.na(country_profile_vars$country_code), ]
countries_and_codes_total <- country_profile_vars[!duplicated(country_profile_vars$country_code) , c("country", "country_code")]internet_usage <- country_profile_vars[ , c("country", "country_code", "Individuals.using.the.Internet..per.100.inhabitants.")]
names(internet_usage) <- c("country", "country_code","internet_usage")
internet_usage[internet_usage$internet_usage == -99, "internet_usage"] <- NA
internet_usage <- correct_names(internet_usage)Echando un vistazo a la variable, salta a la vista que algunos valores están mal introducidos, ya que es imposible que, haciéndolo relastivo a número de personas, en Brasil se acceda más a Internet que en España, pues es bien conocido que es un país con una buena parte de la población viviendo en la pobreza. Por este motivo muy probablemente acabe eliminada.
mobile_cellular_subscriptions <- country_profile_vars[ , c("country" ,"country_code", "Mobile.cellular.subscriptions..per.100.inhabitants.")]
names(mobile_cellular_subscriptions) <- c("country", "country_code","cell_subscr")
mobile_cellular_subscriptions$cell_subscr[mobile_cellular_subscriptions$cell_subscr == -99 | mobile_cellular_subscriptions$cell_subscr == "..."] <- NA
#mobile_cellular_subscriptions[mobile_cellular_subscriptions$country %in% inf_mort_rate$country, ]
mobile_cellular_subscriptions <- correct_names(mobile_cellular_subscriptions)
internet_usage$country <- NULL
mobile_cellular_subscriptions$country <- NULL
mydata <- merge(internet_usage, mobile_cellular_subscriptions, on = "country_code", all = T)gdp <- data.frame(country_profile_vars$country, country_profile_vars$country_code ,country_profile_vars$GDP..Gross.domestic.product..million.current.US..)
names(gdp) <- c("country","country_code", "gdp")
gdp$gdp[gdp$gdp == -99 | gdp$gdp == "..."] <- NA
gdp <- gdp[!is.na(gdp$country), ]
gdp$country <- as.character(gdp$country)
gdp <- correct_names(gdp)
gdp$country <- NULL
mydata <- merge(mydata, gdp, on = "country_code", all = T)library(readr)
worldbankdata <- read_csv("worldbank1.csv")## Parsed with column specification:
## cols(
## Time = col_integer(),
## `Time Code` = col_character(),
## `Country Name` = col_character(),
## `Country Code` = col_character(),
## `GDP per capita (Current US$)` = col_character(),
## `Labor force, female (% of total labor force)` = col_character(),
## `Inflation, consumer prices (annual %)` = col_character(),
## `GNI per capita, Atlas method (current US$)` = col_character(),
## `Maternity leave (days paid)` = col_character(),
## `Ratio of female to male labor force participation rate (%) (modeled ILO estimate)` = col_character(),
## `Contributing family workers, female (% of female employment) (modeled ILO estimate)` = col_character()
## )worldbankdata <- worldbankdata[ , c("Country Name", "GDP per capita (Current US$)", "Labor force, female (% of total labor force)", "Inflation, consumer prices (annual %)", "GNI per capita, Atlas method (current US$)", "Maternity leave (days paid)", "Ratio of female to male labor force participation rate (%) (modeled ILO estimate)", "Contributing family workers, female (% of female employment) (modeled ILO estimate)")]
names(worldbankdata) <- c("country", "gdp_per_cap", "lab_force_fem", "inflation", "gni_per_cap", "matern_leave", "fem_to_male_part_rate", "contr_fam_workers")
worldbankdata[worldbankdata == ".."] <- NA
worldbankdata <- correct_names(worldbankdata)
worldbankdata$country_code <- countrycode(worldbankdata$country, "country.name", "iso3c")## Warning in countrycode(worldbankdata$country, "country.name", "iso3c"): Some values were not matched unambiguously: Caribbean small states, Channel Islands, Eswatini, Euro area, Kosovo, North America, Saint Martin, Small states, South Asia, Worldworldbankdata <- add_code_absent(worldbankdata)
worldbankdata <- worldbankdata[!is.na(worldbankdata$country_code), ]
worldbankdata <- data.frame(worldbankdata)
worldbankdata[ , 2:8] <- apply(worldbankdata[ , 2:8], 2, as.numeric)
worldbankdata$country <- NULL
mydata <- merge(mydata, worldbankdata, on = "country_code", all = T)balance_of_payments <- country_profile_vars[ , c("country", "country_code","Balance.of.payments..current.account..million.US..")]
names(balance_of_payments) <- c("country","country_code", "BoP")
balance_of_payments$BoP[balance_of_payments$BoP == -99 | balance_of_payments$BoP == "..."] <- NA
balance_of_payments <- correct_names(balance_of_payments)
balance_of_payments$country <- NULL
mydata <- merge(mydata, balance_of_payments, on = "country_code", all = T)urban_pop <- country_profile_vars[, c("country", "country_code","Urban.population....of.total.population.")]
names(urban_pop) <- c("country", "country_code","urban_p")
urban_pop$urban_p[urban_pop$urban_p == -99 | urban_pop$urban_p == ".."] <- NA
urban_pop <- correct_names(urban_pop)
urban_pop$country <- NULL
mydata <- merge(mydata, urban_pop, on = "country_code", all = T)refugees <- country_profile_vars[ , c("country", "country_code","Refugees.and.others.of.concern.to.UNHCR..in.thousands." )]
names(refugees) <- c("country", "country_code","refugees")
refugees$refugees[refugees$refugees == -99 | refugees$refugees == ".."] <- NA
refugees <- correct_names(refugees)
refugees[refugees == "~0.0"] <- NA
refugees$country <- NULL
mydata <- merge(mydata,refugees, on = "country_code", all = T)gini_index <- read_csv("gini_index.csv")
names(gini_index)## [1] "Series Name" "Series Code" "Country Name" "Country Code"
## [5] "1990 [YR1990]" "2000 [YR2000]" "2008 [YR2008]" "2009 [YR2009]"
## [9] "2010 [YR2010]" "2011 [YR2011]" "2012 [YR2012]" "2013 [YR2013]"
## [13] "2014 [YR2014]" "2015 [YR2015]" "2016 [YR2016]" "2017 [YR2017]"country_and_code <- gini_index[ , c("Country Name", "Country Code")]
names(country_and_code) <- c("country", "country_code")
gini_index<- gini_index[ ,9:ncol(gini_index)]
gini_index <- apply(gini_index, 2, as.numeric)
gini_index$mean <- rowMeans(gini_index, na.rm = T)
gini_index <- data.frame(cbind(country_and_code, gini_index$mean))
names(gini_index) <- c("country", "country_code", "gini_index")
gini_index <- correct_names(gini_index)
gini_index$country <- NULL
mydata <- merge(mydata, gini_index, on = "country_code", all = T)health <- country_profile_vars[, c("country", "country_code","Health..Total.expenditure....of.GDP.")]
names(health) <- c("country", "country_code","health")
health$health[health$health == -99 | health$health == ".."] <- NA
health <- correct_names(health)
health$country <- NULL
mydata <- merge(mydata,health, on = "country_code", all = T)education <- country_profile_vars[ , c("country", "country_code","Education..Government.expenditure....of.GDP.")]
names(education) <- c("country", "country_code","education")
education$education[education$education == -99 | education$education == ".." | education$education == "..."] <- NA
education <- correct_names(education)
education$country <- NULL
mydata <- merge(mydata,education, on = "country_code", all = T)co2 <- country_profile_vars[ , c("country", "country_code","CO2.emission.estimates..million.tons.tons.per.capita.")]
names(co2) <- c("country", "country_code","co2")
co2$co2[co2$co2 == -99 | co2$co2 == ".."] <- NA
co2 <- correct_names(co2)
co2$country <- NULL
mydata <- merge(mydata,co2, on = "country_code", all = T)productivity <- read_csv("productivity_wb.csv") #as measured by GDP per person employed. ## Parsed with column specification:
## cols(
## Time = col_integer(),
## `Time Code` = col_character(),
## `Country Name` = col_character(),
## `Country Code` = col_character(),
## `GDP per person employed (constant 2011 PPP $)` = col_character()
## )productivity <- productivity[ , c(3, 4, 5)]
names(productivity) <- c("country", "country_code", "productivity")
productivity$productivity[productivity$productivity == -99 | productivity$productivity == ".."] <- NA
productivity <- correct_names(productivity)
productivity <- add_code_absent(productivity)
productivity$country <- NULL
mydata <- merge(mydata, productivity, on = "country_code", all = T)innovation_and_fin_health <- read_csv("innovation_and_debt.csv")## Parsed with column specification:
## cols(
## Time = col_integer(),
## `Time Code` = col_character(),
## `Country Name` = col_character(),
## `Country Code` = col_character(),
## `Scientific and technical journal articles` = col_character(),
## `Patent applications, nonresidents` = col_character(),
## `Patent applications, residents [IP.PAT.RESD]` = col_character(),
## `Trademark applications, total` = col_character(),
## `High-technology exports (current US$)` = col_character(),
## `Charges for the use of intellectual property, receipts (BoP, current US$)` = col_character(),
## `Depth of credit information index (0=low to 8=high)` = col_character(),
## `Depth of the food deficit (kilocalories per person per day)` = col_character(),
## `IMF charges (INT, current US$)` = col_character(),
## `Interest forgiven (current US$)` = col_character()
## )head(innovation_and_fin_health)## # A tibble: 6 x 14
## Time `Time Code` `Country Name` `Country Code` `Scientific and~
## <int> <chr> <chr> <chr> <chr>
## 1 2016 YR2016 Afghanistan AFG 80.4
## 2 2016 YR2016 Albania ALB 191.4
## 3 2016 YR2016 Algeria DZA 4447.1
## 4 2016 YR2016 American Samoa ASM ..
## 5 2016 YR2016 Andorra AND 7.9
## 6 2016 YR2016 Angola AGO 39.1
## # ... with 9 more variables: `Patent applications, nonresidents` <chr>,
## # `Patent applications, residents [IP.PAT.RESD]` <chr>, `Trademark
## # applications, total` <chr>, `High-technology exports (current
## # US$)` <chr>, `Charges for the use of intellectual property, receipts
## # (BoP, current US$)` <chr>, `Depth of credit information index (0=low
## # to 8=high)` <chr>, `Depth of the food deficit (kilocalories per person
## # per day)` <chr>, `IMF charges (INT, current US$)` <chr>, `Interest
## # forgiven (current US$)` <chr>innovation_and_fin_health <- innovation_and_fin_health[ , 3:ncol(innovation_and_fin_health)]
names(innovation_and_fin_health) <- c("country", "country_code", "sci_jour", "patent_nr", "patent_r", "trademark app", "hightech_exp", "int_prop", "credit_inf", "food_deficit", "imf_charges", "int_forgiven")
innovation_and_fin_health[innovation_and_fin_health == ".." | innovation_and_fin_health == -99] <- NA
innovation_and_fin_health <- correct_names(innovation_and_fin_health)
innovation_and_fin_health <- add_code_absent(innovation_and_fin_health)
innovation_and_fin_health$country <- NULL
mydata <- merge(mydata, innovation_and_fin_health, on = "country_code", all = T)renewable_energy <- read_csv("renewable_energy.csv")## Parsed with column specification:
## cols(
## Time = col_integer(),
## `Time Code` = col_character(),
## `Country Name` = col_character(),
## `Country Code` = col_character(),
## `Renewable energy share of TFEC (%)` = col_character(),
## `Renewable electricity share of total electricity output (%)` = col_character()
## )renewable_energy <- renewable_energy[ , 3:ncol(renewable_energy)]
names(renewable_energy) <- c("country", "country_code", "renew_en", "renew_electr")
renewable_energy <- correct_names(renewable_energy)
renewable_energy <- add_code_absent(renewable_energy)
renewable_energy$country <- NULL
mydata <- merge(mydata, renewable_energy, on = "country_code", all = T)library(WHO)
library(rgho)
library(plyr)
library(dplyr)
require(tidyr)
infant_mortality_rate <- data.frame(get_data("MDG_0000000001"))
inf_mort_rate <- infant_mortality_rate %>%
group_by(country) %>%
group_by(year, add = T) %>%
filter(year >= 2010) %>%
select(country, value) %>%
filter(!is.na(country)) %>%
spread(year, value)
inf_mort_rate$inf_mort_rate <- inf_mort_rate$`2017`
inf_mort_rate <- inf_mort_rate[ , c("country", "inf_mort_rate")]
inf_mort_rate <- correct_names(inf_mort_rate)
inf_mort_rate$country_code <- countrycode(inf_mort_rate$country, "country.name", "iso3c")
inf_mort_rate <- add_code_absent(inf_mort_rate)sanitation_services <- read.csv("sanitation_services.csv")
names(sanitation_services) <- c("country", "sanitation_services")
sanitation_services <- sanitation_services[!is.na(sanitation_services$country), ]
sanitation_services$country <- as.character(sanitation_services$country)
sanitation_services <- correct_names(sanitation_services)
sanitation_services$country_code <- countrycode(sanitation_services$country, "country.name", "iso3c")
sanitation_services <- add_code_absent(sanitation_services)
inf_mort_rate$country <- NULL
sanitation_services$country <- NULL
who_df <- merge(inf_mort_rate, sanitation_services, on = "country_code", all = T)deaths_for_household_pollution <- read.csv("deaths_for_household_pollution.csv", sep = ",", fileEncoding = "utf-8", stringsAsFactors = F)
deaths_for_household_pollution$Value <- as.numeric(gsub(' .*', "", deaths_for_household_pollution$Value))
deaths_for_household_pollution <- deaths_for_household_pollution[deaths_for_household_pollution$Cause == "Total", ]
names(deaths_for_household_pollution) <- c("country", "cause", "deaths_household_pollution")
deaths_for_household_pollution$cause <- NULL
deaths_for_household_pollution <- deaths_for_household_pollution[!is.na(deaths_for_household_pollution$country), ]
deaths_for_household_pollution <- correct_names(deaths_for_household_pollution)
deaths_for_household_pollution$country_code <- countrycode(deaths_for_household_pollution$country, "country.name", "iso3c")
deaths_for_household_pollution <- add_code_absent(deaths_for_household_pollution)
deaths_for_household_pollution$country <- NULL
who_df <- merge(who_df, deaths_for_household_pollution, on = "country_code", all = T)water <- read.csv("water.csv", stringsAsFactors = F)
names(water) <- c("country", "water")
water <- water[!is.na(water$country), ]
water <- correct_names(water)
water$country_code <- countrycode(water$country, "country.name", "iso3c")
water <- add_code_absent(water)
water$country <- NULL
who_df <- merge(who_df, water, on = "country_code", all = T)require(tidyr)
library(plyr)
require(dplyr)
energy_wb <- read_csv("energy_wb.csv")## Parsed with column specification:
## cols(
## Time = col_character(),
## `Time Code` = col_character(),
## `Country Name` = col_character(),
## `Country Code` = col_character(),
## `Access to electricity (% of total population) [1.1_ACCESS.ELECTRICITY.TOT]` = col_character(),
## `Access to Clean Fuels and Technologies for cooking (% of total population) [2.1_ACCESS.CFT.TOT]` = col_character()
## )energy_wb <- energy_wb[1:460, c(1, 3:6)]
names(energy_wb) <- c("year", "country", "country_code", "electricity", "clean_fuels")
library(data.table) ## v >= 1.9.6##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, lastenergy_wb <- dcast(setDT(energy_wb), country ~ year, value.var = c("electricity", "clean_fuels"))
energy_wb[energy_wb == ".." | energy_wb == -99] <- NA
energy_new <- data.frame(energy_wb)
for(column in c("electricity", "clean_fuels")) {
nombre_real <- paste0(column, "_2016")
for (i in 1:nrow(energy_wb)) {
if(is.na(energy_new[i, nombre_real])) {
if(!is.na(energy_new[i, paste0(column, "_2015")])) {
energy_new[i, nombre_real] <- energy_new[i, paste0(column, "_2015")]
} else {
next
}
}
}
}
energy_wb <- correct_names(energy_wb)
energy_new$country_code <- countrycode(energy_new$country, "country.name", "iso3c")## Warning in countrycode(energy_new$country, "country.name", "iso3c"): Some values were not matched unambiguously: Channel Islands, Kosovo, Netherlands Antillesenergy_new <- add_code_absent(energy_new)
energy_new <- energy_new[!is.na(energy_new$country_code), ]
energy_new$electricity_2015 <- NULL
energy_new$clean_fuels_2015 <- NULL
energy_new$country <- NULL
mydata <- merge(mydata, energy_new, on = "country_code", all = T)mydata <- merge(mydata, who_df, on = "country_code", all = T)
#codes_and_countries <- mydata[ , c("country", "country_code")]
#mydata$country_code <- NULLEn este chunk de abajo se irán eliminando algunas variables que, después de haber realizado el análisis de factores completo y de haber visualizado suficientemente la relación multidimensional con el resto de variables, resultan no tener mucho que aportar a la construcción del modelo, y hacen que este no pueda captar la mayoría de la varianza total de la muestra de variables.
additional_un_vars <- country_profile_vars[ , c("country_code", "Fertility.rate..total..live.births.per.woman.", "Life.expectancy.at.birth..females.males..years.", "Urban.population.growth.rate..average.annual...", "International.trade..Balance..million.US..", "Employment..Agriculture....of.employed.", "Economy..Agriculture....of.GVA.", "Economy..Industry....of.GVA.")]
additional_un_vars[additional_un_vars== ".." | additional_un_vars == "..." | additional_un_vars == -99 | additional_un_vars == ".../..."] <- NA
names(additional_un_vars) <- c("country_code", "fert_rate", "life_exp", "urban_pop_growth", "int_tr_bal", "emp_agr", "eco_agr", "eco_ind")
life_exp <- c()
for (i in 1:nrow(additional_un_vars)) {
if(!is.na(additional_un_vars[i, "life_exp"])) {
fem <- as.numeric(gsub("/.*", "", additional_un_vars[i, "life_exp"]))
mas <- as.numeric(gsub(".*/", "", additional_un_vars[i, "life_exp"]))
media <- (fem + mas) / 2
life_exp <- c(life_exp, media)
} else {
life_exp <- c(life_exp, NA)
}
}
additional_un_vars$life_exp <- life_exp
mydata <- merge(mydata, additional_un_vars, on = "country_code", all = T)vulner_empl <- read_csv("vulnerable_employment.csv")## Parsed with column specification:
## cols(
## Time = col_character(),
## `Time Code` = col_character(),
## `Country Name` = col_character(),
## `Country Code` = col_character(),
## `Vulnerable employment, total (% of total employment) (modeled ILO estimate) [SL.EMP.VULN.ZS]` = col_character()
## )vulner_empl <- vulner_empl[1:217 , 3:ncol(vulner_empl)]
names(vulner_empl) <- c("country", "country_code", "vulner_empl")
vulner_empl[vulner_empl == ".." | vulner_empl == -99] <- NA
vulner_empl$country <- NULL
mydata <- merge(mydata, vulner_empl, on = "country_code", all = T)La variable domestic credit viene de World Bank y es una variable que nos ayuda a entender el nivel de confianza entre los agentes económicos de un país. Esperamos que cuanto mayor sea el nivel de spread, o cuanto mejor posicionados estén los países en la escala de spread, mayor será también la confianza, identificada en este caso como el credito doméstico al sector privado como porcentaje del GDP en el año 2017.
dom_credit <- read_csv("dom_credit.csv")## Parsed with column specification:
## cols(
## Time = col_integer(),
## `Time Code` = col_character(),
## `Country Name` = col_character(),
## `Country Code` = col_character(),
## `Domestic credit to private sector (% of GDP) [FS.AST.PRVT.GD.ZS]` = col_character()
## )dom_credit <- dom_credit[ , 3:5]
names(dom_credit) <- c("country", "country_code", "dom_credit")
dom_credit[dom_credit == "..." | dom_credit == ".." | dom_credit == -99] <- NA
dom_credit$country <- NULL
mydata <- merge(mydata, dom_credit, on = "country_code", all = T)La variable account ownership es una variable que nos indica el porcentaje de personas de un país que tienen acceso, bien sea por tener una cuenta en una entidad financiera o bien por gozar de este servicio vía móvil, al sistema económico y financiero; cuántas personas forman parte de ese sistema dentro de un país. Intuitivamente, en los países que mejor estén aprovechando este progreso tecnológico habrá actividad económica, y esta variable nos dice cuánta gente participa en esa actividad económica; este es un asunto comentado por los autores de especial importancia a la hora de ir “beyond gdp”, pues para que el bounty sea sostenible, y puedas posicionarte como país arriba en la escala de spread, este bounty tiene que estar bien distribuido entre tu poblacón; de otra forma no existirá equilibrio social ni económico.
acc_own <- read_csv("account_ownership.csv")## Parsed with column specification:
## cols(
## Time = col_integer(),
## `Time Code` = col_character(),
## `Country Name` = col_character(),
## `Country Code` = col_character(),
## `Account ownership at a financial institution or with a mobile-money-service provider (% of population ages 15+) [FX.OWN.TOTL.ZS]` = col_character()
## )acc_own <- acc_own[ , 3:5]
names(acc_own) <- c("country", "country_code", "acc_own")
acc_own[acc_own == "..." | acc_own == ".." | acc_own == -99] <- NA
acc_own$country <- NULL
mydata <- merge(mydata, acc_own, on = "country_code", all = T)educ_staff_comp <- read_csv("educ_staff_comp.csv")## Parsed with column specification:
## cols(
## `Series Name` = col_character(),
## `Series Code` = col_character(),
## `Country Name` = col_character(),
## `Country Code` = col_character(),
## `2010 [YR2010]` = col_character(),
## `2011 [YR2011]` = col_character(),
## `2012 [YR2012]` = col_character(),
## `2013 [YR2013]` = col_character(),
## `2014 [YR2014]` = col_character(),
## `2015 [YR2015]` = col_character(),
## `2016 [YR2016]` = col_character(),
## `2017 [YR2017]` = col_character()
## )educ_staff_comp <- educ_staff_comp[ , 3:ncol(educ_staff_comp)]
educ_staff_comp$ed_comp <- rowMeans(apply(educ_staff_comp[ ,3:ncol(educ_staff_comp)], 2, as.numeric ), na.rm=T)## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción
## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción
## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción
## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción
## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción
## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción
## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coerción
## Warning in apply(educ_staff_comp[, 3:ncol(educ_staff_comp)], 2,
## as.numeric): NAs introducidos por coercióneduc_staff_comp <- educ_staff_comp[ , c(1,2,11)]
names(educ_staff_comp) <- c("country", "country_code", "ed_comp")
educ_staff_comp[educ_staff_comp == "..." | educ_staff_comp == ".." | educ_staff_comp == -99] <- NA
educ_staff_comp <- educ_staff_comp[!is.na(educ_staff_comp$country), ]
educ_staff_comp$country <- NULL
mydata <- merge(mydata, educ_staff_comp, on = "country_code", all = T)prop_miss <- function(x){sum(is.na(x))/length(x)*100}
apply(mydata,2,prop_miss)## country_code internet_usage
## 0.000000 3.361345
## cell_subscr gdp
## 11.344538 11.764706
## gdp_per_cap lab_force_fem
## 21.008403 21.848739
## inflation gni_per_cap
## 31.092437 22.689076
## matern_leave fem_to_male_part_rate
## 25.210084 21.428571
## contr_fam_workers BoP
## 21.428571 23.949580
## urban_p refugees
## 3.361345 31.092437
## gini_index health
## 42.857143 19.327731
## education co2
## 36.974790 11.344538
## productivity sci_jour
## 22.268908 16.806723
## patent_nr patent_r
## 51.680672 54.621849
## trademark app hightech_exp
## 42.857143 43.697479
## int_prop credit_inf
## 49.579832 19.747899
## food_deficit imf_charges
## 50.840336 47.478992
## int_forgiven renew_en
## 47.478992 2.521008
## renew_electr electricity_2016
## 2.521008 10.084034
## clean_fuels_2016 inf_mort_rate
## 20.588235 17.647059
## sanitation_services deaths_household_pollution
## 18.487395 23.529412
## water fert_rate
## 18.487395 10.504202
## life_exp urban_pop_growth
## 10.924370 3.361345
## int_tr_bal emp_agr
## 10.924370 18.067227
## eco_agr eco_ind
## 12.605042 11.764706
## vulner_empl dom_credit
## 20.588235 41.596639
## acc_own ed_comp
## 39.495798 44.537815apply(mydata,1,prop_miss)## [1] 54.166667 14.583333 16.666667 70.833333 2.083333 47.916667 93.750000
## [8] 18.750000 2.083333 2.083333 83.333333 37.500000 8.333333 6.250000
## [15] 4.166667 14.583333 6.250000 6.250000 87.500000 8.333333 2.083333
## [22] 2.083333 12.500000 18.750000 6.250000 2.083333 12.500000 56.250000
## [29] 0.000000 0.000000 16.666667 18.750000 14.583333 6.250000 31.250000
## [36] 8.333333 8.333333 89.583333 6.250000 6.250000 41.666667 8.333333
## [43] 16.666667 20.833333 58.333333 2.083333 18.750000 12.500000 2.083333
## [50] 31.250000 85.416667 70.833333 6.250000 6.250000 6.250000 16.666667
## [57] 39.583333 6.250000 6.250000 4.166667 4.166667 6.250000 33.333333
## [64] 81.250000 6.250000 8.333333 10.416667 6.250000 10.416667 81.250000
## [71] 6.250000 75.000000 43.750000 14.583333 6.250000 2.083333 2.083333
## [78] 79.166667 12.500000 79.166667 12.500000 10.416667 29.166667 10.416667
## [85] 29.166667 62.500000 2.083333 81.250000 68.750000 10.416667 14.583333
## [92] 31.250000 2.083333 6.250000 14.583333 6.250000 2.083333 83.333333
## [99] 4.166667 6.250000 8.333333 18.750000 8.333333 8.333333 6.250000
## [106] 8.333333 4.166667 10.416667 0.000000 2.083333 0.000000 6.250000
## [113] 39.583333 45.833333 4.166667 4.166667 18.750000 10.416667 12.500000
## [120] 12.500000 31.250000 18.750000 62.500000 4.166667 12.500000 6.250000
## [127] 8.333333 6.250000 43.750000 91.666667 4.166667 56.250000 2.083333
## [134] 0.000000 14.583333 0.000000 47.916667 10.416667 10.416667 6.250000
## [141] 8.333333 10.416667 2.083333 79.166667 4.166667 12.500000 72.916667
## [148] 79.166667 4.166667 6.250000 2.083333 81.250000 12.500000 52.083333
## [155] 6.250000 14.583333 12.500000 77.083333 8.333333 8.333333 4.166667
## [162] 60.416667 14.583333 16.666667 2.083333 4.166667 2.083333 6.250000
## [169] 45.833333 20.833333 8.333333 47.916667 43.750000 8.333333 8.333333
## [176] 35.416667 52.083333 16.666667 93.750000 2.083333 2.083333 2.083333
## [183] 14.583333 18.750000 8.333333 12.500000 81.250000 14.583333 4.166667
## [190] 0.000000 43.750000 33.333333 83.333333 4.166667 27.083333 12.500000
## [197] 20.833333 6.250000 10.416667 6.250000 10.416667 70.833333 27.083333
## [204] 29.166667 70.833333 14.583333 8.333333 0.000000 12.500000 87.500000
## [211] 22.916667 16.666667 25.000000 14.583333 2.083333 6.250000 52.083333
## [218] 93.750000 10.416667 2.083333 6.250000 10.416667 10.416667 18.750000
## [225] 91.666667 18.750000 25.000000 70.833333 68.750000 4.166667 18.750000
## [232] 85.416667 14.583333 79.166667 18.750000 2.083333 12.500000 4.166667get_noisy_countries <- function(df) {
require(lava)
if ("country" %ni% names(df)) {
noisy_countries <- c()
for(i in 1:nrow(df)) {
if(sum(is.na(df[i, ])) / ncol(df) > 0.3) {
noisy_countries <- c(noisy_countries, df[i, "country_code"])
} else {
next
}
}
} else {
noisy_countries <- c()
for(i in 1:nrow(df)) {
if(sum(is.na(df[i, ])) / ncol(df) > 0.3) {
noisy_countries <- c(noisy_countries, df[i, "country"])
} else {
next
}
}
}
return(noisy_countries)
}
noisy_c <- get_noisy_countries(mydata)## Loading required package: lava## lava version 1.6.3##
## Attaching package: 'lava'## The following object is masked from 'package:dplyr':
##
## varsnoisy_c## [1] "ABW" "AIA" "AND" "ANT" "ASM" "ATG" "BES" "BMU" "CAF" "CHI" "CIV"
## [12] "COK" "CUB" "CUW" "CYM" "DMA" "ERI" "ESH" "FLK" "FRO" "FSM" "GIB"
## [23] "GLP" "GRL" "GUF" "GUM" "HKG" "IMN" "KIR" "KNA" "LBY" "LIE" "MAC"
## [34] "MAF" "MCO" "MHL" "MNP" "MSR" "MTQ" "MYT" "NCL" "NIU" "NRU" "PLW"
## [45] "PRI" "PRK" "PSE" "PYF" "REU" "SHN" "SMR" "SOM" "SPM" "SXM" "TCA"
## [56] "TKL" "TUV" "TWN" "VAT" "VGB" "VIR" "WLF" "XKX"Vamos a probar eliminando primero las variables que tengan muchos valores ausentes, y después volvemos a ver cuántos países ruidosos tenemos.
remove_noisy_vars <- function(df) {
cols_pos <- c()
for(i in 1:ncol(df)) {
if(sum(is.na(df[, i])) / nrow(df) > 0.4) {
cols_pos <- c(cols_pos, i)
} else {
next
}
}
return(cols_pos)
}
cols_pos <- remove_noisy_vars(mydata)
mydata <- mydata[, -cols_pos]
noisy_c <- get_noisy_countries(mydata)
noisy_c## [1] "ABW" "AIA" "AND" "ANT" "ASM" "BES" "BMU" "CHI" "CIV" "COK" "CUW"
## [12] "CYM" "DMA" "ESH" "FLK" "FRO" "GIB" "GLP" "GRL" "GUF" "GUM" "IMN"
## [23] "KNA" "LIE" "MAC" "MAF" "MCO" "MHL" "MNP" "MSR" "MTQ" "MYT" "NCL"
## [34] "NIU" "NRU" "PLW" "PRI" "PRK" "PYF" "REU" "SHN" "SMR" "SPM" "SXM"
## [45] "TCA" "TKL" "TUV" "TWN" "VAT" "VGB" "VIR" "WLF" "XKX"library(VIM)
aggr_plot <- aggr(mydata, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(mydata), cex.axis=.7, gap=3, ylab=c("Histogram of missing data","Pattern"))
##
## Variables sorted by number of missings:
## Variable Count
## acc_own 0.39495798
## education 0.36974790
## inflation 0.31092437
## refugees 0.31092437
## matern_leave 0.25210084
## BoP 0.23949580
## deaths_household_pollution 0.23529412
## gni_per_cap 0.22689076
## productivity 0.22268908
## lab_force_fem 0.21848739
## fem_to_male_part_rate 0.21428571
## contr_fam_workers 0.21428571
## gdp_per_cap 0.21008403
## clean_fuels_2016 0.20588235
## vulner_empl 0.20588235
## credit_inf 0.19747899
## health 0.19327731
## sanitation_services 0.18487395
## water 0.18487395
## emp_agr 0.18067227
## inf_mort_rate 0.17647059
## sci_jour 0.16806723
## eco_agr 0.12605042
## gdp 0.11764706
## eco_ind 0.11764706
## cell_subscr 0.11344538
## co2 0.11344538
## life_exp 0.10924370
## int_tr_bal 0.10924370
## fert_rate 0.10504202
## electricity_2016 0.10084034
## internet_usage 0.03361345
## urban_p 0.03361345
## urban_pop_growth 0.03361345
## renew_en 0.02521008
## renew_electr 0.02521008
## country_code 0.00000000mydata_df <- mydata
#mydata_numeric <- data.frame(sapply(mydata[ , -c(1,2)], as.numeric))
mydata_df <- mydata_df[!duplicated(mydata_df$country_code), ]
mydata_df[ , 2:ncol(mydata_df)] <- sapply(mydata_df[ , 2:ncol(mydata_df)], as.numeric)Ahora vamos a ver cómo quedaría el missing data plot usando sólo aquellos países con menos de 30% de valores ausentes.
require(lava)
mydata2 <- mydata_df[mydata_df$country_code %ni% noisy_c, ]
aggr_plot <- aggr(mydata2, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(mydata), cex.axis=.7, gap=3, ylab=c("Histogram of missing data","Pattern"))
##
## Variables sorted by number of missings:
## Variable Count
## education 0.256830601
## acc_own 0.229508197
## refugees 0.153005464
## inflation 0.147540984
## BoP 0.092896175
## matern_leave 0.065573770
## productivity 0.049180328
## gni_per_cap 0.043715847
## gdp_per_cap 0.038251366
## lab_force_fem 0.038251366
## fem_to_male_part_rate 0.032786885
## contr_fam_workers 0.032786885
## clean_fuels_2016 0.027322404
## deaths_household_pollution 0.027322404
## emp_agr 0.027322404
## vulner_empl 0.027322404
## health 0.021857923
## sanitation_services 0.016393443
## water 0.016393443
## co2 0.010928962
## sci_jour 0.010928962
## credit_inf 0.010928962
## inf_mort_rate 0.010928962
## urban_pop_growth 0.010928962
## eco_agr 0.005464481
## country_code 0.000000000
## internet_usage 0.000000000
## cell_subscr 0.000000000
## gdp 0.000000000
## urban_p 0.000000000
## renew_en 0.000000000
## renew_electr 0.000000000
## electricity_2016 0.000000000
## fert_rate 0.000000000
## life_exp 0.000000000
## int_tr_bal 0.000000000
## eco_ind 0.000000000Eliminamos algunas variables que hemos visto que, o no contienen información relevante para el análisis (hemos visto al hacer la matriz de correlaciones que no tienen demasiada correlación con el resto de variables) o hemos tenido que descartar por contener demasiados valores ausentes.
mydata2$internet_usage <- NULL
mydata2$matern_leave <- NULL
mydata2$gdp <- NULL
mydata2$BoP <- NULL
#mydata2$cell_subscr <- NULL
mydata2$contr_fam_workers <- NULL
mydata2$refugees <- NULL
#mydata2$health <- NULL
mydata2$education <- NULL
#mydata2$co2 <- NULL
#mydata2$renew_electr <- NULL
mydata2$int_forgiven <- NULL
mydata2$imf_charges <- NULL
mydata2$int_prop <- NULL
mydata2$food_deficit <- NULL
mydata2$patent_nr <- NULL
mydata2$patent_r <- NULL
#mydata2$sci_jour <- NULL
mydata2$trademark.app <- NULL
mydata2$hightech_exp <- NULL
mydata2$int_tr_bal <- NULL
mydata2$eco_ind <- NULL
mydata2$fem_to_male_part_rate <- NULL
mydata2$lab_force_fem <- NULL
mydata2$co2 <- NULL
mydata2$health <- NULL
mydata2$sci_jour <- NULL
mydata2$cell_subscr <- NULLVemos que mejora considerablemente, por lo que probablemente merezca la pena quedarnos solo con 180 observaciones.
if(!require(mice)) {
install.packages("mice")
}## Loading required package: mice## Loading required package: lattice##
## Attaching package: 'mice'## The following object is masked from 'package:tidyr':
##
## complete## The following objects are masked from 'package:base':
##
## cbind, rbindcountries_and_regions <- country_profile_vars[ , c("country_code","country", "Region")]
names(countries_and_regions) <- c("country_code", "country", "region")
mydata2 <- merge(mydata2, countries_and_regions, on = "country_code", all.x = T)
mydata2$region <- as.factor(mydata2$region)
mydata2$country_code <- NULL
countries <- mydata2$country
mydata2$country <- NULL
mydata_imputed <- mice(mydata2, m=1, method='cart', printFlag=FALSE)## Warning: Number of logged events: 35mydata_complete <- mice::complete(mydata_imputed)
sum(sapply(mydata_complete, function(x) { sum(is.na(x)) }))## [1] 0Las variables con las que nos quedamos son las que vienen a continuación. No explicamos todas porque la mayoría son autoexplicativas con el nombre, y por no ser demasiado redundantes, pues el significado de las variables se comenta de forma tangencial a lo largo de todo el análisis.
- Gdp Per Capita: la sacamos de WB
- Inflation: la sacamos de WB. Medida como índice de precios al consumo.
- Gni Per Capita: la sacamos de WB
- Urban Population: porcentaje de población urbana. La sacamos de UN.
- Employment Agriculture: porcentaje de empleo en agricultura. UN
- Economy Agriculture: Importancia de la agricultura en la economía. UN.
- Urban Population Growth: de 2016 a 2017. UN.
- Credit Information: Del 1 al 8 la información de crédito que tiene la gente de un país. WB.
- Productivity: GDP/Hour worked. WB.
- Renewable Energy: Uso de energía renovable. WB.
- Renewable Electricity: Porcentaje de la energía eléctrica que es renovable. WB.
- Electricity 2016: WB.
- Clean Fuels 2016: Acceso a “clean fuels”. WB.
- Infant Mortality Rate: WHO.
- Sanitation Services: WHO.
- Deaths for Household Pollution: WHO.
- Water: WHO.
- Fertility Rate: UN.
- Account Ownership: WB.
- Life Expectancy: UN.
- Vulnerable Employment: WB.
2. Análisis Descriptivo
2.1. Individual
if(!require(fitdistrplus)){
install.packages("fitdistrplus")
}## Loading required package: fitdistrplus## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: survival## Loading required package: npsurv## Loading required package: lseiregions <- mydata_complete$region
mydata_complete$region <- NULL #ya no lo necesitamos
datos_brutos <- data.frame(apply(mydata_complete, 2, scale))
for(i in 1:ncol(datos_brutos)) {
require(ggplot2)
if(is.numeric(datos_brutos[ , i])) {
hist(datos_brutos[ , i], main = names(datos_brutos)[i], col = "blue")
boxplot(datos_brutos[ , i], col = "green", main = names(datos_brutos)[i])
if(is.integer(datos_brutos[ , i])) {
descdist(datos_brutos[ , i], discrete = T) #si todos los elementos del vector son enteros, es una variable discreta que toma suficientes valores.
} else {
descdist(datos_brutos[ , i], discrete = F) #sino, es una variable continua
}
} else if(is.factor(datos_brutos[ , i])) {
counts <- table(datos_brutos[ , i])
barplot(counts, main = names(datos_brutos)[i])
descdist(as.numeric(datos_brutos[ , i]), discrete = T)
} else {
next
}
}## Loading required package: ggplot2##
## Attaching package: 'ggplot2'## The following object is masked from 'package:lava':
##
## vars
Para no alargar demasiado unas explicaciones que tampoco tienen una relevancia crucial en todos los casos en la consecución del análisis multivariante, en esta sección nos centraremos principalmente en visualizar la distribución de todas las variables (escaladas), de tal forma que nos podremos hacer una idea de los outliers que tienen (podrían sesgar el análisis posterior). Para ello, dibujaremos el boxplot e histograma de cada una de las variables, así como, ayudándonos de la librería fitdistrplus, la distribución a la que más se asemejan las mismas. Esto nos permite, de forma sencilla y visual, la curtosis y asimetría que tienen cada una de nuestras variables, y poder compararlas respecto a las esperables para cada una de las distribuciones que incluye (ver leyenda de los gráficos). El gdp per capita se distribuye más o menos como una beta; podemos ver que tiene una fuerte asimetría a la derecha, y, viendo el boxplot, hay varias observaciones que serían outliers en una normal. La variable inflación presenta rarezas, ya que casi todos los países están en torno al 0, mientras que vemos incluso una observación en 12.31 (1231% de inflación). Esto muestra que la mayoría de países tienen una inflación tremendamente similar, siendo muy pequeño el IQR, mientras que algunos países tienen especiales problemas de inflación y contienen outliers. Estas observaciones pueden resultar problemáticas a la hora de aplicar una técnica de reducción de la dimensión como Principal Component Analysis (PCA) y Factor Analysis (FA). Debido a la enorme magnitud de esa observación, vamos a eliminarla. Lo comentado para el gdp per capita parece cumplirse también para otras variables como el gni per capita, productivity, eco_agr (esta variable representa la importancia de la agricultura como sector económico dentro de la economía doméstica). Hay otras variables que también parecen comportarse como una beta, aunque con distintos grados tanto de curtosis como de asimetría. Algo que observamos en muchas variables es que parece haber diferentes “montañas” dentro de la distribución, mostrando que puede haber distintos grupos dentro de los países, dentro de los cuales los datos sí fluctúan más o menos al rededor de una media pero que están separadas por “valles”, es decir países intermedios, siendo el número de estos bastante escaso. Un ejemplo de eso es la variable credit_inf, es decir información sobre crédito de la población. Podemos ver que hay dos medias, una para el mundo desarrollado y otra para los países en vías de desarrollo. Hay muchos países en el que el nivel de conocimiento de servicios crediticios etc es muy inferior al que podemos tener en determinados países de Europa o en América del Norte. Muy relacionada con esta variable está acc_own, es decir el porcentaje de la población que tiene cuenta bancaria o similar (explicación más profunda más arriba). Y también en esta variable podemos ver estas dos “montañas” que antes comentaba, esa separación de la muestra en dos grupos claramente diferenciados. Fert_rate es otra variable que tiene un comportamiento de este tipo, aunque mucho menos pronunciado. Las variables que se comportan de esta forma son, a priori, las que mejor nos van a ayudar a distinguir el nivel de spread, o el posicionamiento con respecto al evento spread (es decir qué tal posicionado está el país con respecto a los demás en términos de calidad de vida derivada del progreso económico y tecnológico). El resto de variables siguen distribuciones o bien similares a la uniforme, o bien similares a la beta con distintos grados de curtosis y también de asimetría (además, tenemos asimetrías positivas y negativas dentro de la muestra, no todas las variables la tienen en la misma dirección). Aparte de los ya comentados outliers de inflation, no observamos ningún otro desajuste excesivo en los datos, por lo que podemos proceder a continuar con la parte más importante del análisis descriptivo, que es la multivariante. Más que cómo están distribuidas nuestras variables, nos interesa saber cómo se distribuyen unas con respecto a las otras; queremos ver si se adivinan correlaciones lineales entre las variables. En definitiva el objetivo es tener un criterio previo que nos defina qué variables pueden o pueden no ser adecuadas y significativas para el análisis posterior, con el fin de eliminar ruido de la muestra y tener una idea previa de la cantidad de información común entre las variables. Cuanta más relación contengan, sin llegar a tener información redundante (necesitamos algo de multicolinearidad para que se puedan extraer factores significativos, pero ésta no puede ser perfecta porque no podríamos utilizar el método necesario para la obtención de factores – matriz no invertible).
datos_brutos <- datos_brutos[!datos_brutos$inflation == max(datos_brutos$inflation), ]2.2. Multivariante
require(GGally)
#datos_brutos <- sapply(mydata[ , -1], as.numeric)
#datos_brutos$inflation <- NULL
#datos_brutos$BoP <- NULL
cormatrix <- cor(datos_brutos)
ggcorr(cormatrix, label = T) 
ggpairs(datos_brutos)
En los gráficos superiores se pueden apreciar la matriz de correlación de nuestras variables. En primer lugar, tenemos este gráfico de la matriz de correlación simple, en el que, como podemos ver, las variables tienen una correlación tremendamente alta unas con otras (tanto positiva como negativa). Aunque en muchos casos pone 1, es porque lo aproxima, y en el siguiente gráfico podemos ver que no llegan nunca a ser 1. Nos hemos asegurado antes de continuar de que no teníamos ningún caso de multicolinearidad perfecta, por lo comentado en el apartado anterior. En el segundo de estos gráficos, además, tenemos el scatterplot de unas variables con respecto a otras, y en la diagonal la pdf, en la que podemos ver cómo se distribuye cada variable, similar a los histogramas realizados antes. Podemos ver que muchas variables tienen relación lineal con las demás.
boxplot(datos_brutos, col="blue", las=2)
library(igraph)##
## Attaching package: 'igraph'## The following objects are masked from 'package:lava':
##
## compare, path## The following object is masked from 'package:tidyr':
##
## crossing## The following objects are masked from 'package:dplyr':
##
## as_data_frame, groups, union## The following objects are masked from 'package:stats':
##
## decompose, spectrum## The following object is masked from 'package:base':
##
## unionR <- cor(datos_brutos)
# visualize size of correlations
heatmap(R, symm = TRUE)
# better to reorder colors
heatmap(sqrt(2*(1-R)), symm = TRUE) # using a distance through correlation
# Correlations induce some network structure
graph.f<-graph.adjacency(R, weighted=TRUE, mode="undirected", diag=F)
plot.igraph(graph.f, edge.width=2,vertex.size=20)
# remove correlations lower than than percentile 20%
g=delete.edges(graph.f, E(graph.f)[ abs(weight) < quantile(abs(R), probs=0.2)])
plot.igraph(g, edge.width=2,vertex.size=20)
# But the relevant multivatiate information is in the inverse of R
iR <- solve(R)
igraph.f<-graph.adjacency(iR, weighted=TRUE, mode="lower", diag=F)
# remove partial correlations lower than perentile 20%
ig <- delete.edges(igraph.f, E(igraph.f)[ abs(weight) < quantile(abs(iR), probs=0.2) ])
plot.igraph(ig, edge.width=2,vertex.size=20)
En este chunk de arriba vemos diferentes formas de visualizar la información que contiene nuestro set de datos, en concreto la información relativa a cómo se relacionan unas variables con otras.
2.3. Inferencia
En este apartado realizaremos algunos test de hipótesis; queremos ver si nuestros datos siguen una distribución normal, entre otras cosas.
datos_brutos = as.matrix(datos_brutos)
n = dim(datos_brutos)[1]
p = dim(datos_brutos)[2]
n## [1] 182p## [1] 21# Check whether mu=0 (vector hypothesis)
mu0 = rep(0,1,p)
Sigma0 <- (1/n) * t(datos_brutos) %*% datos_brutos - mu0 %*% t(mu0)
Sigma0## gdp_per_cap inflation gni_per_cap urban_p
## gdp_per_cap 0.99743869 -0.10779893 0.98791782 0.59331094
## inflation -0.10779893 0.14310958 -0.10909547 -0.08746353
## gni_per_cap 0.98791782 -0.10909547 0.99793090 0.60364108
## urban_p 0.59331094 -0.08746353 0.60364108 0.98543922
## productivity 0.88209569 -0.11348328 0.87177861 0.66279367
## credit_inf 0.25899746 -0.05042613 0.28207298 0.41305309
## renew_en -0.31505805 0.11550050 -0.31936175 -0.48033537
## renew_electr -0.02461287 0.09524137 -0.02939236 -0.17337863
## electricity_2016 0.38913598 -0.12384488 0.39711980 0.53094788
## clean_fuels_2016 0.52316039 -0.09242545 0.53158870 0.64111485
## inf_mort_rate -0.52237889 0.13830966 -0.53320946 -0.55622284
## sanitation_services 0.50727895 -0.10420559 0.51792597 0.57458375
## deaths_household_pollution -0.54700681 0.08785447 -0.55532721 -0.61253540
## water 0.46378713 -0.13528488 0.47426166 0.59812316
## fert_rate -0.47746874 0.12291346 -0.48791983 -0.52964720
## life_exp 0.63628922 -0.14348935 0.64967537 0.62261185
## urban_pop_growth -0.30348990 0.07348143 -0.31300274 -0.30968505
## emp_agr -0.57247374 0.10621856 -0.58193154 -0.69701655
## eco_agr -0.51314458 0.10637209 -0.52011440 -0.59107134
## vulner_empl -0.60612344 0.10856046 -0.61768529 -0.64416840
## acc_own 0.70775770 -0.12814266 0.71954776 0.52866905
## productivity credit_inf renew_en
## gdp_per_cap 0.8820957 0.25899746 -0.3150581
## inflation -0.1134833 -0.05042613 0.1155005
## gni_per_cap 0.8717786 0.28207298 -0.3193618
## urban_p 0.6627937 0.41305309 -0.4803354
## productivity 0.9947228 0.33485030 -0.5035063
## credit_inf 0.3348503 0.98731418 -0.3088054
## renew_en -0.5035063 -0.30880539 0.9996920
## renew_electr -0.1844359 0.06198460 0.5586055
## electricity_2016 0.5193091 0.48992052 -0.7219390
## clean_fuels_2016 0.6426168 0.49230892 -0.7289692
## inf_mort_rate -0.6015568 -0.53903195 0.6212961
## sanitation_services 0.6225482 0.51325062 -0.7252856
## deaths_household_pollution -0.6312588 -0.46800456 0.6701490
## water 0.5673666 0.47216510 -0.6821997
## fert_rate -0.5604317 -0.54391887 0.6201322
## life_exp 0.6750055 0.49890395 -0.6317957
## urban_pop_growth -0.3571955 -0.32136303 0.5096508
## emp_agr -0.6966648 -0.40254688 0.7205614
## eco_agr -0.6337788 -0.45775220 0.6663125
## vulner_empl -0.7287212 -0.43582621 0.6937656
## acc_own 0.7259118 0.45017963 -0.4540027
## renew_electr electricity_2016 clean_fuels_2016
## gdp_per_cap -0.02461287 0.3891360 0.52316039
## inflation 0.09524137 -0.1238449 -0.09242545
## gni_per_cap -0.02939236 0.3971198 0.53158870
## urban_p -0.17337863 0.5309479 0.64111485
## productivity -0.18443594 0.5193091 0.64261682
## credit_inf 0.06198460 0.4899205 0.49230892
## renew_en 0.55860548 -0.7219390 -0.72896921
## renew_electr 0.99446396 -0.2009990 -0.21720152
## electricity_2016 -0.20099898 0.9638105 0.83557168
## clean_fuels_2016 -0.21720152 0.8355717 0.98421657
## inf_mort_rate 0.18683577 -0.7766352 -0.78362738
## sanitation_services -0.19072477 0.8398051 0.86584425
## deaths_household_pollution 0.18878617 -0.7020436 -0.87447132
## water -0.22503654 0.8393339 0.78350552
## fert_rate 0.14273061 -0.7847359 -0.76916617
## life_exp -0.17072089 0.7870422 0.77999648
## urban_pop_growth 0.14962094 -0.5945882 -0.61503664
## emp_agr 0.32390290 -0.7422659 -0.82466812
## eco_agr 0.19220274 -0.6784794 -0.74405649
## vulner_empl 0.26598501 -0.7307746 -0.84512504
## acc_own -0.13166022 0.5272506 0.60586083
## inf_mort_rate sanitation_services
## gdp_per_cap -0.5223789 0.5072789
## inflation 0.1383097 -0.1042056
## gni_per_cap -0.5332095 0.5179260
## urban_p -0.5562228 0.5745838
## productivity -0.6015568 0.6225482
## credit_inf -0.5390319 0.5132506
## renew_en 0.6212961 -0.7252856
## renew_electr 0.1868358 -0.1907248
## electricity_2016 -0.7766352 0.8398051
## clean_fuels_2016 -0.7836274 0.8658443
## inf_mort_rate 0.9777406 -0.8253493
## sanitation_services -0.8253493 0.9749646
## deaths_household_pollution 0.7803206 -0.7905230
## water -0.7905601 0.8450177
## fert_rate 0.8309591 -0.8153358
## life_exp -0.9041731 0.8323763
## urban_pop_growth 0.5805201 -0.6309816
## emp_agr 0.7374044 -0.8024338
## eco_agr 0.6977185 -0.7190240
## vulner_empl 0.7712448 -0.8219014
## acc_own -0.6912545 0.6195550
## deaths_household_pollution water
## gdp_per_cap -0.54700681 0.4637871
## inflation 0.08785447 -0.1352849
## gni_per_cap -0.55532721 0.4742617
## urban_p -0.61253540 0.5981232
## productivity -0.63125881 0.5673666
## credit_inf -0.46800456 0.4721651
## renew_en 0.67014901 -0.6821997
## renew_electr 0.18878617 -0.2250365
## electricity_2016 -0.70204365 0.8393339
## clean_fuels_2016 -0.87447132 0.7835055
## inf_mort_rate 0.78032055 -0.7905601
## sanitation_services -0.79052297 0.8450177
## deaths_household_pollution 0.98807701 -0.6957409
## water -0.69574090 0.9779317
## fert_rate 0.67233025 -0.8176361
## life_exp -0.77115633 0.7972351
## urban_pop_growth 0.48726028 -0.6298389
## emp_agr 0.76058708 -0.7688848
## eco_agr 0.70599806 -0.6559515
## vulner_empl 0.79936960 -0.7424831
## acc_own -0.60697084 0.6063777
## fert_rate life_exp urban_pop_growth
## gdp_per_cap -0.4774687 0.6362892 -0.30348990
## inflation 0.1229135 -0.1434894 0.07348143
## gni_per_cap -0.4879198 0.6496754 -0.31300274
## urban_p -0.5296472 0.6226118 -0.30968505
## productivity -0.5604317 0.6750055 -0.35719553
## credit_inf -0.5439189 0.4989040 -0.32136303
## renew_en 0.6201322 -0.6317957 0.50965075
## renew_electr 0.1427306 -0.1707209 0.14962094
## electricity_2016 -0.7847359 0.7870422 -0.59458823
## clean_fuels_2016 -0.7691662 0.7799965 -0.61503664
## inf_mort_rate 0.8309591 -0.9041731 0.58052005
## sanitation_services -0.8153358 0.8323763 -0.63098164
## deaths_household_pollution 0.6723303 -0.7711563 0.48726028
## water -0.8176361 0.7972351 -0.62983885
## fert_rate 0.9860985 -0.8225073 0.69205254
## life_exp -0.8225073 0.9792918 -0.55331598
## urban_pop_growth 0.6920525 -0.5533160 0.98308739
## emp_agr 0.7156742 -0.7493805 0.60936854
## eco_agr 0.6637815 -0.6713800 0.50745172
## vulner_empl 0.7432691 -0.7550921 0.62703060
## acc_own -0.6857094 0.6957487 -0.49989456
## emp_agr eco_agr vulner_empl acc_own
## gdp_per_cap -0.5724737 -0.5131446 -0.6061234 0.7077577
## inflation 0.1062186 0.1063721 0.1085605 -0.1281427
## gni_per_cap -0.5819315 -0.5201144 -0.6176853 0.7195478
## urban_p -0.6970165 -0.5910713 -0.6441684 0.5286691
## productivity -0.6966648 -0.6337788 -0.7287212 0.7259118
## credit_inf -0.4025469 -0.4577522 -0.4358262 0.4501796
## renew_en 0.7205614 0.6663125 0.6937656 -0.4540027
## renew_electr 0.3239029 0.1922027 0.2659850 -0.1316602
## electricity_2016 -0.7422659 -0.6784794 -0.7307746 0.5272506
## clean_fuels_2016 -0.8246681 -0.7440565 -0.8451250 0.6058608
## inf_mort_rate 0.7374044 0.6977185 0.7712448 -0.6912545
## sanitation_services -0.8024338 -0.7190240 -0.8219014 0.6195550
## deaths_household_pollution 0.7605871 0.7059981 0.7993696 -0.6069708
## water -0.7688848 -0.6559515 -0.7424831 0.6063777
## fert_rate 0.7156742 0.6637815 0.7432691 -0.6857094
## life_exp -0.7493805 -0.6713800 -0.7550921 0.6957487
## urban_pop_growth 0.6093685 0.5074517 0.6270306 -0.4998946
## emp_agr 0.9765394 0.7921508 0.8885272 -0.6508890
## eco_agr 0.7921508 0.9978747 0.7816005 -0.6176567
## vulner_empl 0.8885272 0.7816005 0.9986436 -0.7075213
## acc_own -0.6508890 -0.6176567 -0.7075213 0.9793095S = cov(datos_brutos)
S## gdp_per_cap inflation gni_per_cap urban_p
## gdp_per_cap 1.00293525 -0.10813568 0.99336321 0.59655516
## inflation -0.10813568 0.13916604 -0.10946558 -0.08732963
## gni_per_cap 0.99336321 -0.10946558 1.00343290 0.60694578
## urban_p 0.59655516 -0.08732963 0.60694578 0.99080319
## productivity 0.88694883 -0.11373873 0.87657681 0.66640708
## credit_inf 0.26039690 -0.05012870 0.28360309 0.41526006
## renew_en -0.31679380 0.11604887 -0.32112177 -0.48297745
## renew_electr -0.02476966 0.09614809 -0.02957345 -0.17438613
## electricity_2016 0.39123271 -0.12355619 0.39926602 0.53375447
## clean_fuels_2016 0.52601565 -0.09229357 0.53449408 0.64457316
## inf_mort_rate -0.52522325 0.13831077 -0.53611787 -0.55919643
## sanitation_services 0.51003735 -0.10397210 0.52074767 0.57765276
## deaths_household_pollution -0.54999842 0.08778141 -0.55836787 -0.61584678
## water 0.46630796 -0.13527256 0.47684455 0.60132867
## fert_rate -0.48007372 0.12298955 -0.49058589 -0.53249483
## life_exp 0.63976439 -0.14354615 0.65322858 0.62595576
## urban_pop_growth -0.30513028 0.07322231 -0.31469936 -0.31130931
## emp_agr -0.57559375 0.10602206 -0.58510814 -0.70076535
## eco_agr -0.51599253 0.10719555 -0.52299955 -0.59436766
## vulner_empl -0.60946189 0.10897189 -0.62108866 -0.64770279
## acc_own 0.71162774 -0.12811498 0.72348701 0.53149398
## productivity credit_inf renew_en
## gdp_per_cap 0.8869488 0.26039690 -0.3167938
## inflation -0.1137387 -0.05012870 0.1160489
## gni_per_cap 0.8765768 0.28360309 -0.3211218
## urban_p 0.6664071 0.41526006 -0.4829775
## productivity 1.0001894 0.33665510 -0.5062811
## credit_inf 0.3366551 0.99269887 -0.3105006
## renew_en -0.5062811 -0.31050058 1.0052134
## renew_electr -0.1854848 0.06228075 0.5616989
## electricity_2016 0.5221019 0.49250889 -0.7259092
## clean_fuels_2016 0.6461168 0.49495068 -0.7329845
## inf_mort_rate -0.6048204 -0.54191718 0.6247142
## sanitation_services 0.6259242 0.51598780 -0.7292774
## deaths_household_pollution -0.6347026 -0.47052228 0.6738409
## water 0.5704416 0.47468130 -0.6859544
## fert_rate -0.5634807 -0.54685058 0.6235469
## life_exp 0.6786770 0.50157078 -0.6352723
## urban_pop_growth -0.3591168 -0.32305759 0.5124539
## emp_agr -0.7004523 -0.40467558 0.7245275
## eco_agr -0.6372989 -0.46030991 0.6699983
## vulner_empl -0.7327325 -0.43821117 0.6975950
## acc_own 0.7298646 0.45257730 -0.4564970
## renew_electr electricity_2016 clean_fuels_2016
## gdp_per_cap -0.02476966 0.3912327 0.52601565
## inflation 0.09614809 -0.1235562 -0.09229357
## gni_per_cap -0.02957345 0.3992660 0.53449408
## urban_p -0.17438613 0.5337545 0.64457316
## productivity -0.18548479 0.5221019 0.64611677
## credit_inf 0.06228075 0.4925089 0.49495068
## renew_en 0.56169892 -0.7259092 -0.73298449
## renew_electr 0.99992765 -0.2021877 -0.21845317
## electricity_2016 -0.20218767 0.9689354 0.84005605
## clean_fuels_2016 -0.21845317 0.8400561 0.98956703
## inf_mort_rate 0.18792934 -0.7807692 -0.78785326
## sanitation_services -0.19184354 0.8442786 0.87051810
## deaths_household_pollution 0.18987407 -0.7058076 -0.87922686
## water -0.22634090 0.8438150 0.78773117
## fert_rate 0.14356765 -0.7889475 -0.77333387
## life_exp -0.17172326 0.7912392 0.78420597
## urban_pop_growth 0.15050103 -0.5977366 -0.61834437
## emp_agr 0.32575539 -0.7462058 -0.82911799
## eco_agr 0.19324569 -0.6822764 -0.74819930
## vulner_empl 0.26746968 -0.7347733 -0.84976867
## acc_own -0.13244676 0.5300124 0.60910828
## inf_mort_rate sanitation_services
## gdp_per_cap -0.5252232 0.5100374
## inflation 0.1383108 -0.1039721
## gni_per_cap -0.5361179 0.5207477
## urban_p -0.5591964 0.5776528
## productivity -0.6048204 0.6259242
## credit_inf -0.5419172 0.5159878
## renew_en 0.6247142 -0.7292774
## renew_electr 0.1879293 -0.1918435
## electricity_2016 -0.7807692 0.8442786
## clean_fuels_2016 -0.7878533 0.8705181
## inf_mort_rate 0.9830195 -0.8297788
## sanitation_services -0.8297788 0.9802128
## deaths_household_pollution 0.7845417 -0.7947950
## water -0.7948054 0.8495564
## fert_rate 0.8354528 -0.8197373
## life_exp -0.9090499 0.8368492
## urban_pop_growth 0.5836201 -0.6343540
## emp_agr 0.7413522 -0.8067332
## eco_agr 0.7016113 -0.7230368
## vulner_empl 0.7754755 -0.8264101
## acc_own -0.6949551 0.6228522
## deaths_household_pollution water
## gdp_per_cap -0.54999842 0.4663080
## inflation 0.08778141 -0.1352726
## gni_per_cap -0.55836787 0.4768446
## urban_p -0.61584678 0.6013287
## productivity -0.63470260 0.5704416
## credit_inf -0.47052228 0.4746813
## renew_en 0.67384090 -0.6859544
## renew_electr 0.18987407 -0.2263409
## electricity_2016 -0.70580758 0.8438150
## clean_fuels_2016 -0.87922686 0.7877312
## inf_mort_rate 0.78454171 -0.7948054
## sanitation_services -0.79479505 0.8495564
## deaths_household_pollution 0.99347013 -0.6994952
## water -0.69949516 0.9832127
## fert_rate 0.67597366 -0.8220567
## life_exp -0.77533005 0.8015216
## urban_pop_growth 0.48987387 -0.6332119
## emp_agr 0.76469682 -0.7730071
## eco_agr 0.70992641 -0.6596134
## vulner_empl 0.80376379 -0.7465550
## acc_own -0.61023749 0.6096098
## fert_rate life_exp urban_pop_growth
## gdp_per_cap -0.4800737 0.6397644 -0.30513028
## inflation 0.1229895 -0.1435461 0.07322231
## gni_per_cap -0.4905859 0.6532286 -0.31469936
## urban_p -0.5324948 0.6259558 -0.31130931
## productivity -0.5634807 0.6786770 -0.35911679
## credit_inf -0.5468506 0.5015708 -0.32305759
## renew_en 0.6235469 -0.6352723 0.51245389
## renew_electr 0.1435676 -0.1717233 0.15050103
## electricity_2016 -0.7889475 0.7912392 -0.59773657
## clean_fuels_2016 -0.7733339 0.7842060 -0.61834437
## inf_mort_rate 0.8354528 -0.9090499 0.58362015
## sanitation_services -0.8197373 0.8368492 -0.63435404
## deaths_household_pollution 0.6759737 -0.7753300 0.48987387
## water -0.8220567 0.8015216 -0.63321189
## fert_rate 0.9914698 -0.8269578 0.69579132
## life_exp -0.8269578 0.9845879 -0.55626958
## urban_pop_growth 0.6957913 -0.5562696 0.98842537
## emp_agr 0.7195284 -0.7533990 0.61262517
## eco_agr 0.6674789 -0.6751259 0.51028844
## vulner_empl 0.7473516 -0.7592346 0.63046839
## acc_own -0.6894042 0.6994783 -0.50255305
## emp_agr eco_agr vulner_empl acc_own
## gdp_per_cap -0.5755938 -0.5159925 -0.6094619 0.7116277
## inflation 0.1060221 0.1071956 0.1089719 -0.1281150
## gni_per_cap -0.5851081 -0.5229995 -0.6210887 0.7234870
## urban_p -0.7007654 -0.5943677 -0.6477028 0.5314940
## productivity -0.7004523 -0.6372989 -0.7327325 0.7298646
## credit_inf -0.4046756 -0.4603099 -0.4382112 0.4525773
## renew_en 0.7245275 0.6699983 0.6975950 -0.4564970
## renew_electr 0.3257554 0.1932457 0.2674697 -0.1324468
## electricity_2016 -0.7462058 -0.6822764 -0.7347733 0.5300124
## clean_fuels_2016 -0.8291180 -0.7481993 -0.8497687 0.6091083
## inf_mort_rate 0.7413522 0.7016113 0.7754755 -0.6949551
## sanitation_services -0.8067332 -0.7230368 -0.8264101 0.6228522
## deaths_household_pollution 0.7646968 0.7099264 0.8037638 -0.6102375
## water -0.7730071 -0.6596134 -0.7465550 0.6096098
## fert_rate 0.7195284 0.6674789 0.7473516 -0.6894042
## life_exp -0.7533990 -0.6751259 -0.7592346 0.6994783
## urban_pop_growth 0.6126252 0.5102884 0.6304684 -0.5025531
## emp_agr 0.9818050 0.7965663 0.8934050 -0.6543633
## eco_agr 0.7965663 1.0033760 0.7859281 -0.6211058
## vulner_empl 0.8934050 0.7859281 1.0041535 -0.7114010
## acc_own -0.6543633 -0.6211058 -0.7114010 0.9846057lambda <- n * log(det(Sigma0)/det(S))
lambda## [1] -12.50472pvalue <- 1 - pchisq(lambda,p)
pvalue## [1] 1# Accept mu=0
# Check whether Sigma=diag (independence under normality)
R = cor(datos_brutos)
lambda <- -n * log(det(R))
lambda## [1] 5248.799pvalue <- 1 - pchisq(lambda,p*(p-1)/2)
pvalue## [1] 0Como podemos ver arriba, después de hacer inferencia podemos decir que aceptamos la H0 de que \(\mu = 0\), mientras que descartamos la hipótesis de que nuestras variables estén incorreladas.
S <- cov(datos_brutos)
muhat = colMeans(datos_brutos)
dM2 <- mahalanobis(datos_brutos,muhat,S) # Mahalanobis distances (square)
dM2## 1 2 3 4 5 6
## 28.994962 68.555443 21.388509 8.404907 7.558137 13.180230
## 7 8 9 10 11 12
## 18.374634 9.696618 7.720667 21.828417 40.022056 7.125166
## 13 14 15 16 17 18
## 18.604817 17.877130 24.750341 6.219358 6.096611 22.091780
## 19 20 21 22 23 24
## 23.295159 13.033114 15.661010 14.348899 8.837935 15.919010
## 25 26 27 28 29 30
## 31.224956 25.371622 15.826720 32.180285 6.444384 39.388383
## 31 32 33 34 35 36
## 6.182840 20.168302 15.728870 51.061244 25.523341 11.542334
## 37 38 39 40 41 42
## 23.726673 5.609281 15.090946 15.153698 4.929121 2.919185
## 43 44 45 46 47 48
## 5.552209 32.563770 11.065248 9.515942 13.065805 9.020102
## 49 50 51 52 53 54
## 43.520519 59.652245 6.078852 7.027677 28.643175 7.920348
## 55 56 57 58 59 60
## 25.020633 5.560876 28.249454 42.183325 6.815805 27.789687
## 61 62 63 64 65 66
## 30.990331 15.348570 21.747149 23.384304 53.329767 6.557230
## 67 68 69 70 71 72
## 13.205011 16.021994 13.215927 10.896540 9.845490 7.955181
## 73 74 75 76 77 78
## 34.236913 4.782653 7.203811 22.786987 29.744670 20.300150
## 79 80 81 82 83 84
## 22.821519 35.986102 12.566905 8.842708 6.403739 18.261907
## 85 86 87 88 89 90
## 8.159629 11.226229 30.453753 13.543354 20.859019 33.518235
## 91 92 93 94 95 96
## 67.035594 24.947342 37.321071 14.566408 56.126314 22.700155
## 97 98 99 100 101 102
## 22.080959 24.545888 41.604811 6.617250 140.925812 8.102901
## 103 104 105 106 107 108
## 14.297767 15.066632 21.627115 21.104485 8.191407 11.423334
## 109 110 111 112 113 114
## 21.665455 9.577928 31.399360 19.964891 23.727328 18.253856
## 115 116 117 118 119 120
## 16.961657 7.941349 28.143204 7.602743 43.726340 30.849777
## 121 122 123 124 125 126
## 46.365370 13.530767 6.101376 30.791715 27.825926 11.319838
## 127 128 129 130 131 132
## 43.361454 33.042003 13.166181 9.251674 13.843996 31.744640
## 133 134 135 136 137 138
## 5.289577 5.845897 13.592292 20.383757 34.692212 14.075633
## 139 140 141 142 143 144
## 10.032621 32.351157 23.336812 21.270008 13.903470 14.511892
## 145 146 147 148 149 151
## 18.214657 29.745228 9.861480 46.161559 10.564650 24.408264
## 152 153 154 155 156 157
## 26.068018 4.663597 8.031725 12.554002 29.258330 5.160675
## 158 159 160 161 162 163
## 16.341235 37.561653 21.422667 14.616840 32.627211 20.300257
## 164 165 166 167 168 169
## 34.837545 22.561278 28.966642 6.720782 10.242489 16.788131
## 170 171 172 173 174 175
## 25.786257 17.064537 13.055307 13.931525 21.598127 16.581537
## 176 177 178 179 180 181
## 18.760829 12.186404 23.239589 42.546579 15.428137 20.548363
## 182 183
## 21.887457 21.878487plot(sort(dM2,index.return=TRUE)$x,pch=19,col=grey.colors(n=n,start=0,end=.75),xlab="",ylab="",main="Mahalanobis distances (square)")
## Q-Q plot for Chi^2 data against true theoretical distribution:
qqplot(dM2, qchisq(ppoints(n), df = p),xlab="quantiles of Mahalanobis distances",
ylab="Quantiles for Chi-square with p df",pch=19,col="darkred",
main = expression("Q-Q plot for" ~~ {chi^2}[p]))
qqline(dM2, distribution = function(q) qchisq(q, df = p),
prob = c(0.1, 0.6), col = "blue")
Aquí vemos que no se distribuye de forma demasiado normal. PCA a priori no necesita de ninguna suposición de normalidad para que la extracción de la varianza máxima en los menores componentes posibles (los primeros son los que mayor porcentaje de varianza explican) sea eficiente. FA utiliza la suposición de normalidad en las uniquesses del modelo, ya que \(\epsilon \sim N(0, \sigma)\). El hecho de que nuestros datos no se distribuyan de forma normal podría afectar al comportamiento de estas uniquesses, haciendo que los test que utilizamos para ver si los factores son suficientes o no dejarían de tener sentido, pues el p-value se calcula bajo la suposición de normalidad de los datos. Por este motivo en la práctica a menudo no se hace demasiado caso a ese p-value a la hora de evaluar un modelo concreto, sino que se emplea más para comparar unos modelos con otros y ver cuáles explican mejor.
3. PCA & Factor Analysis
Como tenemos 21 variables en nuestra muestra, existen \(2^{p} \approx 2^{21} = 2097152\) diferentes relaciones. El primer objetivo es, por lo tanto, ser capaces de reducir la dimensión, limpiando el ruido. En alta dimensión casi todas las observaciones son outliers; cuando aumentamos nuestros datos, especialmente el número de variables, aumenta la información pero el ruido lo hace de una forma más rápida. Por eso con cada dimensión adicional que incluyamos en nuestro set de datos estaremos incluyendo mucho más ruido que información.
Inicialmente usaremos PCA para detectar cuánta varianza común tienen nuestras variables, para después tratar de encontrar y definir las variables mencionadas al principio del estudio, en especial la variable spread. Como ya se comentó antes, el spread es en definitiva cómo se distribuye ese bounty, y es una variable latente más fácil de encontrar y de definir que el bounty, pues para estimar este necesitaríamos las variables en forma de serie temporal y, dado que para muchas variables hemos tenido que apañarnos con datos de años pasados al 2017 o con medias de los últimos años (para subsanar la ausencia de observaciones más recientes en muchos países), esto sería una tarea infactible muy probablemente. Por eso nos centraremos principalmente en encontrar el spread, es decir la creciente desigualdad causada en gran parte por el desarrollo tecnológico cuyos frutos están desigualmente distribuidos.
FA es lo mismo que PCA si le quitas las uniquesses (o si estas son 0):
PCA:
\(\Sigma = V \Lambda V^{T}\)
FA:
\(\Sigma = L L^{T} + \Psi\)
Si \(\Psi = 0\), FA es igual que PCA con \(L = V\Lambda^{1/2}\).
Para poder definir formalmente la variable spread, debemos buscarla en forma de variable latente mediante un Factor Analysis, que nos permita además ver cuánta varianza explica y deja de explicar la variable oculta spread. Sin embargo, realizamos el PCA primero para ver con cuántos factores será buena idea dotar al modelo, cuánta varianza común hay entre nuestras variables y la posibilidad de reducir la dimensión eficientemente.
3.1. Principal Component Analysis
pca <- prcomp(datos_brutos, center = T, scale = T)
screeplot(pca,main="Screeplot",col="blue",type="barplot",pch=19)
pca## Standard deviations (1, .., p=21):
## [1] 3.59763732 1.37042017 1.14657322 0.96253675 0.84945791 0.73966812
## [7] 0.67354053 0.64932915 0.55332286 0.54219114 0.46175454 0.44023830
## [13] 0.37831575 0.37058149 0.35487052 0.32479451 0.30399767 0.27949216
## [19] 0.23136380 0.20584373 0.09165425
##
## Rotation (n x k) = (21 x 21):
## PC1 PC2 PC3 PC4
## gdp_per_cap 0.1924576 0.49691736 -0.07803238 0.003498832
## inflation -0.1079339 -0.01727274 0.33725021 0.797897121
## gni_per_cap 0.1952434 0.48950047 -0.06895215 -0.002237555
## urban_p 0.2022804 0.16050288 -0.05623977 0.213630935
## productivity 0.2199064 0.35303008 -0.14848017 0.095576458
## credit_inf 0.1532914 -0.05586961 0.42122225 -0.115073872
## renew_en -0.2102537 0.27506181 0.29521566 -0.142928811
## renew_electr -0.0717081 0.25033208 0.69820706 -0.090506917
## electricity_2016 0.2393336 -0.22448831 0.08555264 -0.058065324
## clean_fuels_2016 0.2539321 -0.11149131 0.05623341 0.167900170
## inf_mort_rate -0.2492296 0.07778861 -0.10960187 0.159272093
## sanitation_services 0.2549754 -0.13437552 0.09746871 0.034300013
## deaths_household_pollution -0.2394578 0.03309485 -0.04174147 -0.185848424
## water 0.2439229 -0.15506465 0.05377820 -0.107947086
## fert_rate -0.2416519 0.12205732 -0.16985839 0.182212338
## life_exp 0.2538097 0.01119038 0.06984229 -0.146461950
## urban_pop_growth -0.1846428 0.20705657 -0.12273970 0.077605863
## emp_agr -0.2521812 0.04464787 0.10093908 -0.191426824
## eco_agr -0.2289405 0.03035463 0.01017343 -0.133410786
## vulner_empl -0.2533910 0.01013596 0.04523106 -0.150156519
## acc_own 0.2195788 0.21986228 0.01935181 -0.157304329
## PC5 PC6 PC7
## gdp_per_cap -0.114164775 0.005083875 -0.11581999
## inflation -0.172176012 -0.212300371 -0.31813459
## gni_per_cap -0.100490013 -0.016017818 -0.11296041
## urban_p 0.383880655 0.133240270 0.27499225
## productivity -0.002476745 -0.052855788 -0.03559855
## credit_inf 0.569439765 -0.557853795 0.02786649
## renew_en -0.057692902 0.008488227 0.06340036
## renew_electr -0.101203832 0.419789420 0.26115582
## electricity_2016 0.025844979 0.240202817 -0.09373285
## clean_fuels_2016 0.017726890 0.171381240 0.00116888
## inf_mort_rate -0.026813855 -0.029303878 0.27622622
## sanitation_services -0.042070660 0.133868600 -0.13727783
## deaths_household_pollution -0.120129005 -0.197163678 0.12829073
## water -0.039932964 0.141911082 -0.09424735
## fert_rate 0.139291134 0.105954085 0.11815650
## life_exp 0.017329116 0.127982450 -0.34642161
## urban_pop_growth 0.621533207 0.284822827 -0.29564473
## emp_agr 0.041977648 -0.020294054 -0.26754542
## eco_agr -0.092278205 0.012240581 -0.50538452
## vulner_empl 0.095839515 0.034318499 -0.21107176
## acc_own -0.130596943 -0.413851902 -0.02313192
## PC8 PC9 PC10
## gdp_per_cap 0.009647005 0.21046911 -0.074492312
## inflation 0.077825079 -0.01170776 0.213860280
## gni_per_cap 0.019504881 0.18312985 -0.084822640
## urban_p 0.653974934 -0.28237279 -0.005637807
## productivity -0.026488158 0.27230759 0.036581697
## credit_inf -0.010439600 0.26430081 -0.216412234
## renew_en 0.143621115 -0.36181378 -0.067124106
## renew_electr -0.110201741 0.08718815 0.048694568
## electricity_2016 0.070694066 0.32346888 0.139484180
## clean_fuels_2016 -0.090140506 0.03667897 -0.317561885
## inf_mort_rate 0.112352458 0.37416017 -0.105206490
## sanitation_services -0.030303241 0.15255269 -0.002555221
## deaths_household_pollution 0.336026447 0.29702805 0.486051063
## water 0.316868790 0.12657818 0.206200369
## fert_rate -0.148193007 0.06991818 -0.230951167
## life_exp 0.040733732 -0.16815145 0.100861539
## urban_pop_growth -0.241526402 -0.05942265 0.357526995
## emp_agr -0.060756207 0.10328084 -0.052958568
## eco_agr 0.405772198 -0.02144044 -0.460955589
## vulner_empl 0.137070945 0.15146356 0.114733881
## acc_own -0.152099675 -0.33988491 0.237723801
## PC11 PC12 PC13
## gdp_per_cap -0.14801344 0.001008133 0.10938078
## inflation -0.04385091 -0.004475410 0.05500664
## gni_per_cap -0.14868511 0.009584271 0.13184416
## urban_p -0.05335836 -0.272317139 -0.07970450
## productivity 0.07605659 0.105223475 -0.24213800
## credit_inf -0.05331763 0.112913890 0.01080246
## renew_en -0.02379598 0.557890932 0.32596461
## renew_electr 0.12373901 -0.200603595 -0.20786441
## electricity_2016 0.14508275 0.072899056 0.17486561
## clean_fuels_2016 0.15368369 -0.057071427 0.11452099
## inf_mort_rate 0.30529794 -0.132938747 0.07541762
## sanitation_services 0.10723198 0.167175678 -0.12344804
## deaths_household_pollution -0.02666143 0.217827160 -0.23058604
## water 0.26135638 0.117241028 0.48049265
## fert_rate 0.03304768 0.135884645 0.35453714
## life_exp -0.35559877 0.003318500 -0.10672764
## urban_pop_growth 0.27424400 0.137588291 -0.08835846
## emp_agr -0.08244732 -0.368030295 0.05243925
## eco_agr 0.38380045 0.046628532 -0.27886733
## vulner_empl -0.23914122 -0.405963699 0.39493736
## acc_own 0.54247332 -0.308554134 0.15171033
## PC14 PC15 PC16
## gdp_per_cap 0.010711120 0.028195710 0.225007438
## inflation 0.000495221 -0.033811711 -0.004010458
## gni_per_cap 0.048615312 0.006540802 0.321001578
## urban_p -0.025124159 -0.067700225 -0.076870334
## productivity -0.176909370 0.083007258 -0.700500594
## credit_inf 0.086189948 -0.007615905 0.047727022
## renew_en -0.354097808 0.047413695 -0.086186809
## renew_electr 0.196196388 -0.042391119 -0.027516007
## electricity_2016 -0.324903114 -0.567810451 0.028543704
## clean_fuels_2016 -0.451291339 -0.068645203 0.056849669
## inf_mort_rate -0.093309449 0.133294401 0.283603157
## sanitation_services 0.210023728 0.200773886 -0.060755524
## deaths_household_pollution -0.028975365 -0.238341946 0.039196148
## water 0.291182873 0.406418265 -0.121903272
## fert_rate 0.484057979 -0.466722652 -0.299024950
## life_exp 0.173032312 -0.248754800 0.060906438
## urban_pop_growth -0.067082657 0.103747810 0.188137640
## emp_agr -0.242219239 0.131110940 -0.311181149
## eco_agr 0.119556888 -0.119734612 0.042662069
## vulner_empl -0.055317974 0.007330463 -0.085152380
## acc_own 0.042899877 -0.230669674 -0.015777446
## PC17 PC18 PC19
## gdp_per_cap 0.002153638 0.090949157 0.097326245
## inflation 0.034938828 0.049119963 -0.007421070
## gni_per_cap -0.044412343 0.139996813 0.069551130
## urban_p -0.150470506 0.105819691 -0.023794930
## productivity 0.184404861 -0.239139516 -0.100347108
## credit_inf 0.055718134 -0.008838201 0.003339975
## renew_en -0.143676048 -0.076769373 -0.171241215
## renew_electr 0.111091441 -0.033903849 0.077475221
## electricity_2016 0.098332618 0.118834906 -0.269289191
## clean_fuels_2016 -0.145724401 -0.090989711 0.520256080
## inf_mort_rate -0.089926692 -0.133066435 -0.528564451
## sanitation_services -0.796782837 -0.147643684 -0.093193205
## deaths_household_pollution -0.205124823 0.072134384 0.371511061
## water 0.246544325 0.175067389 0.092675216
## fert_rate -0.146530578 0.100934166 0.085406315
## life_exp -0.005764918 -0.111388624 -0.344949381
## urban_pop_growth 0.018623575 -0.003567013 0.062840838
## emp_agr -0.254572925 0.625711795 -0.103143970
## eco_agr 0.118648953 -0.076920538 0.067160273
## vulner_empl -0.078661483 -0.610078893 0.125719345
## acc_own -0.142927505 -0.064926494 -0.024837813
## PC20 PC21
## gdp_per_cap -0.138153531 -0.7150331379
## inflation -0.001232173 -0.0006412151
## gni_per_cap -0.068321328 0.6907793980
## urban_p -0.085943661 -0.0078154078
## productivity 0.066614199 0.0741340121
## credit_inf 0.024821085 -0.0197353848
## renew_en 0.007798176 -0.0069109998
## renew_electr 0.016436726 0.0122815014
## electricity_2016 -0.321254707 0.0125601908
## clean_fuels_2016 0.447103471 -0.0071226079
## inf_mort_rate 0.350833901 -0.0340143911
## sanitation_services -0.187958094 -0.0070090088
## deaths_household_pollution 0.228454492 -0.0114428416
## water 0.185898660 -0.0045036233
## fert_rate 0.117816471 -0.0151802643
## life_exp 0.604475818 -0.0461375741
## urban_pop_growth 0.022260596 -0.0032954164
## emp_agr 0.109466182 -0.0053926851
## eco_agr -0.069884202 0.0032495515
## vulner_empl -0.132770875 0.0337197537
## acc_own 0.030706564 -0.0171334248summary(pca)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 3.5976 1.37042 1.1466 0.96254 0.84946 0.73967
## Proportion of Variance 0.6163 0.08943 0.0626 0.04412 0.03436 0.02605
## Cumulative Proportion 0.6163 0.70576 0.7684 0.81248 0.84684 0.87290
## PC7 PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.6735 0.64933 0.55332 0.5422 0.46175 0.44024
## Proportion of Variance 0.0216 0.02008 0.01458 0.0140 0.01015 0.00923
## Cumulative Proportion 0.8945 0.91458 0.92916 0.9432 0.95331 0.96254
## PC13 PC14 PC15 PC16 PC17 PC18
## Standard deviation 0.37832 0.37058 0.3549 0.32479 0.3040 0.27949
## Proportion of Variance 0.00682 0.00654 0.0060 0.00502 0.0044 0.00372
## Cumulative Proportion 0.96935 0.97589 0.9819 0.98691 0.9913 0.99503
## PC19 PC20 PC21
## Standard deviation 0.23136 0.20584 0.09165
## Proportion of Variance 0.00255 0.00202 0.00040
## Cumulative Proportion 0.99758 0.99960 1.00000En primer lugar, utilizamos PCA para ver cómo se comportan en alta dimensión nuestras variables. Vemos que el primer componente es muy fuerte, y consigue captar la mayor parte de la varianza total de la muestra, en concreto 3.59 es la desviación estándar de nuestro primer componente. La del segundo es 1.39, y la del tercero 1.13.Los que vienen después captan cada vez menos varianza. Aunque el primero tiene una gran diferencia con respecto a los otros dos (proporción acumulada de varianza en este primer componente es del 61.35%), podría ser buena idea coger 3, ascendiendo hasta el 76.5% de proporción de varianza acumulada. Esto sugiere que hay muchas relaciones ocultas entre nuestras variables y que éstas se pueden explicar en gran medida mediante el uso de 1 componente, aunque con 3 podemos explicar casi un 80% de la varianza total de la muestra.
dfplot <- data.frame(country = countries[-150], region = regions[-150], pc1 = pca$x[ , 1], pc2 = pca$x[ , 2])
qplot(dfplot$pc1,dfplot$pc2,col=dfplot$region,label=dfplot$country,size=I(.1)) +
labs(title="First two principal components", x="PC1", y="PC2",color="Region") +
theme_bw() + theme(legend.position="bottom") + geom_text(size=2, hjust=.4, vjust=-.4)
En el gráfico superior podemos ver el primer contra el segundo componente en los ejes x e y, con los países y sus regiones que definen el color de la etiqueta. Vemos que claramente países como Luxemburgo, Noruega, Suiza, Islandia o Singapur tienen mucha fuerza en ambos componentes, y son los países que intuitivamente tienen un mayor nivel de vida y por lo tanto esperamos que estén en ese grupo de bien posicionados con respecto a la Second Machine Age y por tanto se estén comiendo una mayor parte de la tarta de los beneficios que el progreso tecnológico provee (bounty).
Como veremos a continuación, los países que aparecen en la parte más alta (alto componente 2) son los más potentes a nivel económico, mientras que los que aparecen más a la derecha alto componente 1) son los mejores en general, los que mejor posicionados están respecto a la escala de spread.
Vamos a dibujar los dos primeros componentes.
barplot(pca$rotation[,1], las=2, col="darkgreen", main = "pc1")
barplot(pca$rotation[,2], las=2, col="darkgreen", main = "pc2")
barplot(pca$rotation[ , 3], las = 2, col = "darkgreen", main = "pc3")
El primer componente vemos que utiliza casi todas las variables, consigue explicar una parte de todas. Tiene peso positivo para el gdp per capita, el gni per capita, urban population, productivity, credit information, electricity 2016, clean fuels 2016, sanitation services, water, life expectancy y account ownership. Por otro lado, tiene pesos negativos para inflation, renewable energy, renewable electricity, infant mortality rate, deaths for household pollution, fertility rate, urban population growth, emp agr (empleo en agricultura – porcentaje del labour force dedicada a la agricultura), economy agriculture (porcentaje de la economía doméstica proveniente de la agricultura) y vulnerable employment. Por lo tanto podríamos decir que este primer componente es algo parecido a lo que venimos buscando, pues explica bastante bien las diferencias entre los países, desde económicas hasta sanitarias pasando por el tipo de energía que utilizan.
El segundo componente tiene un peso alto positivo para gdp per capita, gni per capita, urban population, productivity, energías renovables, electricidad renovable, negativo para electricity 2016, negativo para clean fuels 2016, negativo para el agua (muy bajo; no le da mucho peso ni al acceso al agua ni a los servicios sanitarios), positivo para urban pop growth y positivo para account ownership. Vemos por tanto que este segundo componente es principalmente un componente económico puramente. En él España por ejemplo tiene menos fuerza que otros países como Irlanda (aunque siendo un paraíso fiscal tiene trampa esta conclusión).
Por último, el gráfico 3 nos muestra que el tercer componente contiene principalmente las variables credit information, renewable energy y renewable electricity, es decir, cómo de avanzado está el país en cuanto a sostenibilidad y el nivel de conocimiento general que hay del crédito, lo cual es fundamental para que un país aproveche la ola del progreso tecnológico, pues sin una participación activa del sistema económico, un menor porcentaje de la población tendrá acceso a los beneficios de ese desarrollo. El uso de energías limpias y la cultura general de la población son fundamentales para que el crecimiento económico y ese race against the machine del país sea sostenible.
Después de este PCA, vamos a utilizar FA para finalizar con el estudio y definir formalmente la variable spread, la cual posiblemente encontremos contenida entre varios factores distintos. Vamos a incluir 3 factores, pues como hemos visto arriba el uso de 3 componentes principales nos permitiría explicar el 76% de la varianza de la muestra y, lo que es mejor, los 3 son interpretables fácilmente y es muy sencillo aplicarlos a la cuestión que tratamos de resolver. Viendo el screeplot bien podríamos haber seguido el criterio del codo y haber cortado en 1, pero va a ser insuficiente para explicar las diferencias. Esperamos que el primer factor se encargue de explicar el spread de forma general, atendiendo a todas las variables que esperamos sean explicadas por esta variable latente; los otros dos factores posiblemente, basándonos en los resultados de PCA (si utilizamos FA sin rotación) sean el componente del spread que tiene más que ver con lo puramente económico y el componente del mismo que tiene más que ver con cuestiones de sostenibilidad y participación ciudadana del sistema.
3.2. Factor Analysis
No aplicaremos ninguna rotación, pues tras haber probado con varios tipos nos encontramos con que rotation = “none” es el método que nos da unos factores más interpretables, además de acumular la mayor parte de la varianza captada por los factores en el primer factor, lo cual es positivo para nuestro objetivo principal que es encontrar la variable latente spread.
factors <- factanal(datos_brutos, factors = 3, rotation = "none", scores = "regression", lower = 0.1)
factors##
## Call:
## factanal(x = datos_brutos, factors = 3, scores = "regression", rotation = "none", lower = 0.1)
##
## Uniquenesses:
## gdp_per_cap inflation
## 0.100 0.846
## gni_per_cap urban_p
## 0.100 0.443
## productivity credit_inf
## 0.123 0.664
## renew_en renew_electr
## 0.316 0.811
## electricity_2016 clean_fuels_2016
## 0.164 0.130
## inf_mort_rate sanitation_services
## 0.111 0.114
## deaths_household_pollution water
## 0.252 0.188
## fert_rate life_exp
## 0.177 0.100
## urban_pop_growth emp_agr
## 0.526 0.114
## eco_agr vulner_empl
## 0.311 0.129
## acc_own
## 0.320
##
## Loadings:
## Factor1 Factor2 Factor3
## gdp_per_cap 0.736 0.635
## inflation -0.367 -0.138
## gni_per_cap 0.745 0.622
## urban_p 0.719 0.131 -0.152
## productivity 0.816 0.444 -0.116
## credit_inf 0.520 -0.169 0.192
## renew_en -0.719 0.318 0.257
## renew_electr -0.224 0.173 0.330
## electricity_2016 0.839 -0.361
## clean_fuels_2016 0.899 -0.212 -0.130
## inf_mort_rate -0.889 0.198 -0.242
## sanitation_services 0.907 -0.250
## deaths_household_pollution -0.853 0.117
## water 0.857 -0.267
## fert_rate -0.850 0.239 -0.208
## life_exp 0.917 0.254
## urban_pop_growth -0.634 0.269
## emp_agr -0.895 0.101 0.272
## eco_agr -0.801 0.197
## vulner_empl -0.903 0.225
## acc_own 0.788 0.219 0.101
##
## Factor1 Factor2 Factor3
## SS loadings 12.683 1.732 0.643
## Proportion Var 0.604 0.082 0.031
## Cumulative Var 0.604 0.686 0.717
##
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 768.77 on 150 degrees of freedom.
## The p-value is 8.13e-84cbind(factors$loadings, factors$uniquenesses)## Factor1 Factor2 Factor3
## gdp_per_cap 0.7359662 0.63452938 4.155716e-02 0.1000000
## inflation -0.3665357 -0.02282047 -1.383439e-01 0.8459872
## gni_per_cap 0.7445106 0.62211085 4.887288e-02 0.1000000
## urban_p 0.7187899 0.13091917 -1.523659e-01 0.4429849
## productivity 0.8161107 0.44424276 -1.156209e-01 0.1232401
## credit_inf 0.5196299 -0.16925103 1.921919e-01 0.6644133
## renew_en -0.7188026 0.31770524 2.569784e-01 0.3163405
## renew_electr -0.2236911 0.17329489 3.299529e-01 0.8111107
## electricity_2016 0.8392060 -0.36107547 4.217055e-02 0.1635819
## clean_fuels_2016 0.8989310 -0.21237754 -1.295020e-01 0.1300528
## inf_mort_rate -0.8894618 0.19823458 -2.421372e-01 0.1109326
## sanitation_services 0.9071072 -0.24988358 1.530662e-02 0.1144710
## deaths_household_pollution -0.8526162 0.11689647 8.460251e-02 0.2522163
## water 0.8569504 -0.26721429 7.670451e-02 0.1883512
## fert_rate -0.8501949 0.23862068 -2.082418e-01 0.1768658
## life_exp 0.9171641 -0.06917325 2.536867e-01 0.1000000
## urban_pop_growth -0.6335769 0.26899404 -2.472499e-05 0.5262214
## emp_agr -0.8951666 0.10121577 2.724716e-01 0.1141917
## eco_agr -0.8005125 0.09459981 1.968513e-01 0.3114842
## vulner_empl -0.9034610 0.06565532 2.245987e-01 0.1290051
## acc_own 0.7883023 0.21906972 1.007613e-01 0.3204361Vemos que con estos 3 factores podemos explicar el 71.6% de la varianza total de la muestra, siendo la proporción de varianza que explica cada factor la siguiente:
Factor 1: 60.3%
Factor 2: 8.2%
Factor 3: 3.1%
Viendo las uniquesses de cada variable, la mayoría son bastante bajas, lo cual nos indica que estos 3 factores explican la mayoría de esas variables. Sin embargo, hay algunas más altas como las de inflation (0.84), urban population (0.44), credit information (0.657), renewable electricity (0.81), urban population growth (0.52). Para mejorar el modelo, podríamos eliminar estas variables, pues están introduciendo ruido a la hora de encontrar la variable latente spread.
barplot(factors$loadings[,1], las = 2, col = "darkgreen", main = "factor1")
Vemos que el primer factor sería, como esperábamos, uno que tiene en cuenta todas las variables prácticamente. Podemos identificar este primer factor como el spread, la variable principal que estábamos buscando.
barplot(factors$loadings[,2], las = 2, col = "darkgreen", main = "factor1")
barplot(factors$loadings[,3], las = 2, col = "darkgreen", main = "factor1")
El segundo factor sería el componente económico del spread, es decir cuál es la parte económica del “bounty” (beneficios debidos al progreso tecnológico), mientras que el tercer factor sería el componente calidad de vida más allá de lo económico, pues como vemos tiene un peso fuerte para las variables credit information, renewable energy, renewable electricity, negativo para infant mortality rate, fertility rate; y por último positivo para emp agr, life expectancy, eco agr, vulnerable employment y account ownership. Lo que menos nos cuadraría sería el vulnerable employment teniendo una relación positiva con el factor calidad de vida más allá de lo económico (o sostenibilidad como dijimos antes), pues esperaríamos que tuvieran una relación negativa.
dfplot <- data.frame(country = countries[-150], region = regions[-150], factor1 = factors$scores[ , 1], factor2 = factors$scores[ , 2])
qplot(dfplot$factor1,dfplot$factor2,col=dfplot$region,label=dfplot$country,size=I(.1)) +
labs(title="First two factors", x="SPREAD IN GENERAL (FACTOR 1)", y="ECONOMIC PART OF SPREAD (FACTOR 2)",color="Region") +
theme_bw() + theme(legend.position="bottom") + geom_text(size=2, hjust=.4, vjust=-.4)
Es especialmente llamativo el caso de Luxemburgo, que es claramente la que mejor posicionada está en el spread; tanto en el general como en la parte más económica del mismo. Los países que están más a la derecha y más arriba son los que se están repartiendo el pastel del progreso tecnológico principalmente. Vemos que coincide con los resultados que podríamos esoerar de forma intuitiva, con países muy desarrollados marcando el ritmo de esta carrera por el progreso: Singapur, Noruega, Islandia… Concluimos por tanto que hemos encontrado lo que veníamos buscando, un factor que pudiéramos identificar como la variable spread. La hemos encontrado realmente repartida en 3 factores: spread general, spread tecnológico y spread sostenible.
Probaremos, por último, a realizar el mismo análisis eliminando las variables que tengan uniquesses demasiado altas, para tratar de afinar un poco más la puntería. De momento nos quedaríamos con el modelo presentado hasta ahora.
datos_reduced <- data.frame(datos_brutos)
datos_reduced[ , c("inflation", "urban_p", "credit_inf", "renew_electr", "urban_pop_growth")] <- NULL
factors_reduced <- factanal(datos_reduced, factors = 3, rotation = "none", scores = "regression", lower = 0.1)
factors_reduced##
## Call:
## factanal(x = datos_reduced, factors = 3, scores = "regression", rotation = "none", lower = 0.1)
##
## Uniquenesses:
## gdp_per_cap gni_per_cap
## 0.100 0.100
## productivity renew_en
## 0.126 0.337
## electricity_2016 clean_fuels_2016
## 0.162 0.122
## inf_mort_rate sanitation_services
## 0.107 0.112
## deaths_household_pollution water
## 0.244 0.191
## fert_rate life_exp
## 0.188 0.100
## emp_agr eco_agr
## 0.129 0.307
## vulner_empl acc_own
## 0.127 0.322
##
## Loadings:
## Factor1 Factor2 Factor3
## gdp_per_cap 0.738 0.632
## gni_per_cap 0.747 0.620
## productivity 0.817 0.441 -0.114
## renew_en -0.716 0.316 0.224
## electricity_2016 0.839 -0.365
## clean_fuels_2016 0.899 -0.218 -0.149
## inf_mort_rate -0.891 0.202 -0.244
## sanitation_services 0.908 -0.253
## deaths_household_pollution -0.855 0.125
## water 0.854 -0.269
## fert_rate -0.848 0.235 -0.193
## life_exp 0.918 0.263
## emp_agr -0.891 0.102 0.257
## eco_agr -0.799 0.213
## vulner_empl -0.902 0.234
## acc_own 0.788 0.223
##
## Factor1 Factor2 Factor3
## SS loadings 11.302 1.586 0.445
## Proportion Var 0.706 0.099 0.028
## Cumulative Var 0.706 0.806 0.833
##
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 536.32 on 75 degrees of freedom.
## The p-value is 7.69e-71cbind(factors_reduced$loadings, factors_reduced$uniquesses)## Factor1 Factor2 Factor3
## gdp_per_cap 0.7383910 0.63207775 0.0363344297
## gni_per_cap 0.7466168 0.61993000 0.0434230561
## productivity 0.8168544 0.44062459 -0.1138407175
## renew_en -0.7164421 0.31627830 0.2237601169
## electricity_2016 0.8387445 -0.36521045 0.0306045157
## clean_fuels_2016 0.8989023 -0.21755158 -0.1488108709
## inf_mort_rate -0.8906982 0.20172388 -0.2437547718
## sanitation_services 0.9077532 -0.25303007 0.0003041284
## deaths_household_pollution -0.8545754 0.12484377 0.0988855429
## water 0.8544821 -0.26897616 0.0796860173
## fert_rate -0.8481071 0.23528000 -0.1933136076
## life_exp 0.9182668 -0.07396433 0.2625352882
## emp_agr -0.8911070 0.10239945 0.2569312519
## eco_agr -0.7988382 0.09736057 0.2126371243
## vulner_empl -0.9022263 0.06629560 0.2344317865
## acc_own 0.7878089 0.22293001 0.0884406168Vemos que el modelo funciona bastante mejor, pues los 3 factores conjuntamente explican un 83.1% de la varianza de la muestra, mientras que el primer factor explica el 70.5%. Como dijimos antes, el p-value se usa en la práctica para comparar modelos. Si comparamos el p-value en cada modelo, vemos que en este segundo modelo después de eliminar esas variables que estaban introduciendo ruido en la muestra el p-value baja significativamente. Podríamos quedarnos, pues, con este modelo, pues descartamos que los factores que componen el spread (spread general + spread económico + spread sostenible) sean variables latentes que expliquen la inflación y el resto de variables que hemos eliminado.
Dibujamos los gráficos correspondientes a este último modelo:
dfplot <- data.frame(country = countries[-150], region = regions[-150], factor1 = factors_reduced$scores[ , 1], factor2 = factors_reduced$scores[ , 2])
qplot(dfplot$factor1,dfplot$factor2,col=dfplot$region,label=dfplot$country,size=I(.1)) +
labs(title="First two factors", x="SPREAD IN GENERAL (FACTOR 1)", y="ECONOMIC PART OF SPREAD (FACTOR 2)",color="Region") +
theme_bw() + theme(legend.position="bottom") + geom_text(size=2, hjust=.4, vjust=-.4)
Vemos que Chad, Somalia o Sierra Leona son algunos de los países que tienen en el spread general una puntuación muy negativa, sin embargo en el spread económico en concreto no salen tan mal parados, y tienen una puntuación similar a países como Australia.
barplot(factors_reduced$loadings[,1], las = 2, col = "darkgreen", main = "factor1")
barplot(factors_reduced$loadings[,2], las = 2, col = "darkgreen", main = "factor2")
barplot(factors_reduced$loadings[,3], las = 2, col = "darkgreen", main = "factor3")
Vemos que los factores que nos salen, por la matriz L que estamos dibujando vector a vector arriba (del factor 1 al factor 3), son muy similares a los obtenidos cuando incluíamos todas las variables. La mejora viene del lado de la reducción del ruido. Al trabajar con 16 variables en lugar de con 21 podemos extraer más información de la muestra y así definir más claramente la variable spread.
4. Visualización de Resultados
A continuación dibujaremos los resultados obtenidos sobre un mapa del mundo, mostrando el score de cada país en cada uno de los 3 factores, que recordemos son:
Factor 1: Spread General
Factor 2: Spread Económico
Factor 3: Spread Sostenible
Primero dibujamos los gráficos para el primer modelo:
if(!require(rworldmap)) {
install.packages("rworldmap")
}
matched <- joinCountryData2Map(data.frame(country=countries[-150], value=factors$scores
[,1]), joinCode="NAME", nameJoinColumn="country")## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your datamatched2 <- joinCountryData2Map(data.frame(country=countries[-150], value=factors$scores
[,2]), joinCode="NAME", nameJoinColumn="country")## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your datamatched3 <- joinCountryData2Map(data.frame(country=countries[-150], value=factors$scores
[,3]), joinCode="NAME", nameJoinColumn="country")## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your datamapCountryData(matched, nameColumnToPlot="value", mapTitle="General Spread by Factor Analysis",
catMethod = "pretty", colourPalette = "topo")
mapCountryData(matched2, nameColumnToPlot="value", mapTitle="Economic Spread by Factor Analysis",
catMethod = "pretty", colourPalette = "topo")
mapCountryData(matched3, nameColumnToPlot="value", mapTitle="Sustainable Spread by Factor Analysis",
catMethod = "pretty", colourPalette = "topo")
# zoom in Africa
mapCountryData(matched, nameColumnToPlot="value", mapRegion = "Africa", mapTitle="General Spread by Factor Analysis: Zoom in Africa",
catMethod = "pretty", colourPalette = "topo")
#zoom in Europe
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Europe", mapTitle = "General Spread by Factor Analysis: Zoom in Europe", catMethod="pretty", colourPalette = "topo")
#zoom in America
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "North America", mapTitle = "General Spread by Factor Analysis: Zoom in North America", catMethod="pretty", colourPalette = "topo")
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Latin America", mapTitle = "General Spread by Factor Analysis: Zoom in South America", catMethod="pretty", colourPalette = "topo")
#zoom in Asia
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Asia", mapTitle = "General Spread by Factor Analysis: Zoom in Asia", catMethod="pretty", colourPalette = "topo")
Vemos que el factor verdaderamente revelador es el del general spread; es el que mejor nos consigue explicar las diferencias entre países y es por lo tanto el resultado al que más atención deberemos prestar en nuestro estudio. Podemos ver que el mundo occidental en general, con América del Norte, Australia, Europa (Western), y algunas regiones aisladas (Japón, United Arab Emirates), son los países que tienen un mejor posicionamiento con respecto a los cambios que está sufriendo el mundo debido al progreso tecnológico, y por lo tanto son los países que más valor están absorbiendo del que se está creando (bounty). Vemos que en definitiva estas dos variables (spread y bounty) están muy relacionadas; el spread es el que mide cómo se distribuye el bounty en el mundo, es decir es una escala que nos permite determinar qué países están aventajando a los demás debido al progreso, qué países están bien posicionados en la balanza cada vez menos equilibrada de la calidad de vida relativa en el mundo. Vemos que existe una desigualdad estructural demográfica, lo cual nos viene a decir que los países colonialistas del siglo pasado son también, en su mayoría, los que están liderando el cambio tecnológico, por lo tanto distanciándose más del resto de países. Especialmente en África observamos que es un continente que se está quedando atrasado respecto al progreso, ya que muchos países de este contienente tienen un score verdaderamente bajo en el spread general. Este mapa nos sirve también para relacionar esta desigualdad derivada del progreso tecnológico con la situación política internacional. Vemos que algunos países como Israel, pese a estar en una zona en la que la mayoría de los países tienen un score mediano aproximadamente en el factor 1, tienen un score equiparable al de Europa, EEUU o Australia, debido a sus conexiones políticas internacionales.
Los mejores países para vivir en líneas generales serían Irlanda (que se cuela entre los países con un score más alto en spread general), Noruega, Suiza, Luxemburgo, Singapur… Estos son los países que mejor puntúan en el spread general, y por lo tanto los que más se están distanciando del resto.
En cuanto al spread económico puro, podemos ver que destacan muy pocos países en este aspecto. Estos son Suiza, Noruega, Islandia e Irlanda. De forma negativa vemos a Egipto y a otros países del norte de África, y algunos países de Europa del Este.
Respecto al Spread Sostenible, vemos que el país que mejor puntúa sobre el resto de forma notable es Nepal, lo cual encaja debido a su religión; precisamente el budismo se centra en buscar el equilibrio; por lo que tiene sentido que también crezcan de forma sostenible (con pocas muertes infantiles, acceso al agua y a la sanidad, alta esperanza de vida, alto uso de energías renovables…). Además, tener un porcentaje importante de economía agrícola podría verse como una forma de crecer más sostenible, pues es una economía enfocada en el alimento, además de acarrear generalmente una mayor calidad de la comida doméstica y ser ésta más saludable.
Dibujaremos por último los mismos gráficos para el segundo modelo de factores, aunque esperamos que los resultados sean prácticamente los mismos como se explica en la anterior sección, pues cambia muy poco la interpretación de los mismos.
matched <- joinCountryData2Map(data.frame(country=countries[-150], value=factors_reduced$scores
[,1]), joinCode="NAME", nameJoinColumn="country")## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your datamatched2 <- joinCountryData2Map(data.frame(country=countries[-150], value=factors_reduced$scores
[,2]), joinCode="NAME", nameJoinColumn="country")## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your datamatched3 <- joinCountryData2Map(data.frame(country=countries[-150], value=factors_reduced$scores
[,3]), joinCode="NAME", nameJoinColumn="country")## 176 codes from your data successfully matched countries in the map
## 6 codes from your data failed to match with a country code in the map
## 67 codes from the map weren't represented in your datamapCountryData(matched, nameColumnToPlot="value", mapTitle="General Spread by Factor Analysis",
catMethod = "pretty", colourPalette = "topo")## You asked for 7 categories, 10 were used due to pretty() classification
mapCountryData(matched2, nameColumnToPlot="value", mapTitle="Economic Spread by Factor Analysis",
catMethod = "pretty", colourPalette = "topo")
mapCountryData(matched3, nameColumnToPlot="value", mapTitle="Sustainable Spread by Factor Analysis",
catMethod = "pretty", colourPalette = "topo")
# zoom in Africa
mapCountryData(matched, nameColumnToPlot="value", mapRegion = "Africa", mapTitle="General Spread by Factor Analysis: Zoom in Africa",
catMethod = "pretty", colourPalette = "topo")## You asked for 7 categories, 10 were used due to pretty() classification
#zoom in Europe
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Europe", mapTitle = "General Spread by Factor Analysis: Zoom in Europe", catMethod="pretty", colourPalette = "topo")## You asked for 7 categories, 10 were used due to pretty() classification
#zoom in America
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "North America", mapTitle = "General Spread by Factor Analysis: Zoom in North America", catMethod="pretty", colourPalette = "topo")## You asked for 7 categories, 10 were used due to pretty() classification
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Latin America", mapTitle = "General Spread by Factor Analysis: Zoom in South America", catMethod="pretty", colourPalette = "topo")## You asked for 7 categories, 10 were used due to pretty() classification
#zoom in Asia
mapCountryData(matched, nameColumnToPlot = "value", mapRegion = "Asia", mapTitle = "General Spread by Factor Analysis: Zoom in Asia", catMethod="pretty", colourPalette = "topo")## You asked for 7 categories, 10 were used due to pretty() classification